pipe_control.diffusion

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

43 -def bmrb_read(star):

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

140 141

142 -def bmrb_write(star):

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

266 267 268

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

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

306 307

308 -def data_names():

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

315 316

317 -def delete():

318 """Function for deleting diffusion tensor data.""" 319 320 # Test if the current data pipe exists. 321 check_pipe() 322 323 # Test if diffusion tensor data exists. 324 if not diff_data_exists(): 325 raise RelaxNoTensorError('diffusion') 326 327 # Delete the diffusion data. 328 del(cdp.diff_tensor)

329 330

331 -def diff_data_exists(pipe=None):

332 """Function for determining if diffusion data exists in the current data pipe. 333 334 @param pipe: The data pipe to search for data in. 335 @type pipe: str 336 @return: The answer to the question. 337 @rtype: bool 338 """ 339 340 # The data pipe to check. 341 if pipe == None: 342 pipe = pipes.cdp_name() 343 344 # Get the data pipe. 345 dp = pipes.get_pipe(pipe) 346 347 # Test if the data structure exists. 348 if hasattr(dp, 'diff_tensor'): 349 return True 350 else: 351 return False

352 353

354 -def display():

355 """Function for displaying the diffusion tensor.""" 356 357 # Test if the current data pipe exists. 358 check_pipe() 359 360 # Test if diffusion tensor data exists. 361 if not diff_data_exists(): 362 raise RelaxNoTensorError('diffusion') 363 364 # Printout. 365 if cdp.diff_tensor.type == 'sphere': 366 tensor_info_table(type='sphere', tm=cdp.diff_tensor.tm, Diso=cdp.diff_tensor.Diso, fixed=cdp.diff_tensor.fixed) 367 elif cdp.diff_tensor.type == 'spheroid': 368 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.Dper, Dratio=cdp.diff_tensor.Dratio, theta=cdp.diff_tensor.theta, phi=cdp.diff_tensor.phi, fixed=cdp.diff_tensor.fixed) 369 elif cdp.diff_tensor.type == 'ellipsoid': 370 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)

371 372

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

374 """Function for setting up a ellipsoidal diffusion tensor. 375 376 @param params: The diffusion tensor parameter. 377 @type params: float 378 @param time_scale: The correlation time scaling value. 379 @type time_scale: float 380 @param d_scale: The diffusion tensor eigenvalue scaling value. 381 @type d_scale: float 382 @param angle_units: The units for the angle parameters which can be either 'deg' or 'rad'. 383 @type angle_units: str 384 @param param_types: The type of parameters supplied. These correspond to 0: {tm, Da, theta, 385 phi}, 1: {Diso, Da, theta, phi}, 2: {tm, Dratio, theta, phi}, 3: {Dpar, 386 Dper, theta, phi}, 4: {Diso, Dratio, theta, phi}. 387 @type param_types: int 388 """ 389 390 # The diffusion type. 391 cdp.diff_tensor.set_type('ellipsoid') 392 393 # Set the parameter values. 394 if params[0] != None: 395 # (tm, Da, Dr, alpha, beta, gamma). 396 if param_types == 0: 397 # Unpack the tuple. 398 tm, Da, Dr, alpha, beta, gamma = params 399 400 # Scaling. 401 tm = tm * time_scale 402 Da = Da * d_scale 403 404 # Set the parameters. 405 set(value=[tm, Da, Dr], param=['tm', 'Da', 'Dr']) 406 407 # (Diso, Da, Dr, alpha, beta, gamma). 408 elif param_types == 1: 409 # Unpack the tuple. 410 Diso, Da, Dr, alpha, beta, gamma = params 411 412 # Scaling. 413 Diso = Diso * d_scale 414 Da = Da * d_scale 415 416 # Set the parameters. 417 set(value=[Diso, Da, Dr], param=['Diso', 'Da', 'Dr']) 418 419 # (Dx, Dy, Dz, alpha, beta, gamma). 420 elif param_types == 2: 421 # Unpack the tuple. 422 Dx, Dy, Dz, alpha, beta, gamma = params 423 424 # Scaling. 425 Dx = Dx * d_scale 426 Dy = Dy * d_scale 427 Dz = Dz * d_scale 428 429 # Set the parameters. 430 set(value=[Dx, Dy, Dz], param=['Dx', 'Dy', 'Dz']) 431 432 # (Dxx, Dyy, Dzz, Dxy, Dxz, Dyz). 433 elif param_types == 3: 434 # Unpack the tuple. 435 Dxx, Dyy, Dzz, Dxy, Dxz, Dyz = params 436 437 # Build the tensor. 438 tensor = zeros((3, 3), float64) 439 tensor[0, 0] = Dxx 440 tensor[1, 1] = Dyy 441 tensor[2, 2] = Dzz 442 tensor[0, 1] = tensor[1, 0] = Dxy 443 tensor[0, 2] = tensor[2, 0] = Dxz 444 tensor[1, 2] = tensor[2, 1] = Dyz 445 446 # Scaling. 447 tensor = tensor * d_scale 448 449 # Eigenvalues. 450 Di, R, alpha, beta, gamma = tensor_eigen_system(tensor) 451 452 # Set the parameters. 453 set(value=[Di[0], Di[1], Di[2]], param=['Dx', 'Dy', 'Dz']) 454 455 # Change the angular units. 456 angle_units = 'rad' 457 458 # Unknown parameter combination. 459 else: 460 raise RelaxUnknownParamCombError('param_types', param_types) 461 462 # Convert the angles to radians. 463 if angle_units == 'deg': 464 print("Converting the angles from degrees to radian units.") 465 alpha = (alpha / 360.0) * 2.0 * pi 466 beta = (beta / 360.0) * 2.0 * pi 467 gamma = (gamma / 360.0) * 2.0 * pi 468 469 # Set the orientational parameters. 470 set(value=[alpha, beta, gamma], param=['alpha', 'beta', 'gamma'])

471 472

473 -def fold_angles(sim_index=None):

474 """Wrap the Euler or spherical angles and remove the glide reflection and translational symmetries. 475 476 Wrap the angles such that:: 477 478 0 <= theta <= pi, 479 0 <= phi <= 2pi, 480 481 and:: 482 483 0 <= alpha <= 2pi, 484 0 <= beta <= pi, 485 0 <= gamma <= 2pi. 486 487 488 For the simulated values, the angles are wrapped as:: 489 490 theta - pi/2 <= theta_sim <= theta + pi/2 491 phi - pi <= phi_sim <= phi + pi 492 493 and:: 494 495 alpha - pi <= alpha_sim <= alpha + pi 496 beta - pi/2 <= beta_sim <= beta + pi/2 497 gamma - pi <= gamma_sim <= gamma + pi 498 499 500 @param sim_index: The simulation index. If set to None then the actual values will be folded 501 rather than the simulation values. 502 @type sim_index: int or None 503 """ 504 505 506 # Wrap the angles. 507 ################## 508 509 # Spheroid. 510 if cdp.diff_tensor.type == 'spheroid': 511 # Get the current angles. 512 theta = cdp.diff_tensor.theta 513 phi = cdp.diff_tensor.phi 514 515 # Simulated values. 516 if sim_index != None: 517 theta_sim = cdp.diff_tensor.theta_sim[sim_index] 518 phi_sim = cdp.diff_tensor.phi_sim[sim_index] 519 520 # Normal value. 521 if sim_index == None: 522 cdp.diff_tensor.set(param='theta', value=wrap_angles(theta, 0.0, 2.0*pi)) 523 cdp.diff_tensor.set(param='phi', value=wrap_angles(phi, 0.0, 2.0*pi)) 524 525 # Simulated theta and phi values. 526 else: 527 cdp.diff_tensor.set(param='theta', value=wrap_angles(theta_sim, theta - pi, theta + pi), category='sim', sim_index=sim_index) 528 cdp.diff_tensor.set(param='phi', value=wrap_angles(phi_sim, phi - pi, phi + pi), category='sim', sim_index=sim_index) 529 530 # Ellipsoid. 531 elif cdp.diff_tensor.type == 'ellipsoid': 532 # Get the current angles. 533 alpha = cdp.diff_tensor.alpha 534 beta = cdp.diff_tensor.beta 535 gamma = cdp.diff_tensor.gamma 536 537 # Simulated values. 538 if sim_index != None: 539 alpha_sim = cdp.diff_tensor.alpha_sim[sim_index] 540 beta_sim = cdp.diff_tensor.beta_sim[sim_index] 541 gamma_sim = cdp.diff_tensor.gamma_sim[sim_index] 542 543 # Normal value. 544 if sim_index == None: 545 cdp.diff_tensor.set(param='alpha', value=wrap_angles(alpha, 0.0, 2.0*pi)) 546 cdp.diff_tensor.set(param='beta', value= wrap_angles(beta, 0.0, 2.0*pi)) 547 cdp.diff_tensor.set(param='gamma', value=wrap_angles(gamma, 0.0, 2.0*pi)) 548 549 # Simulated alpha, beta, and gamma values. 550 else: 551 cdp.diff_tensor.set(param='alpha', value=wrap_angles(alpha_sim, alpha - pi, alpha + pi), category='sim', sim_index=sim_index) 552 cdp.diff_tensor.set(param='beta', value= wrap_angles(beta_sim, beta - pi, beta + pi), category='sim', sim_index=sim_index) 553 cdp.diff_tensor.set(param='gamma', value=wrap_angles(gamma_sim, gamma - pi, gamma + pi), category='sim', sim_index=sim_index) 554 555 556 # Remove the glide reflection and translational symmetries. 557 ########################################################### 558 559 # Spheroid. 560 if cdp.diff_tensor.type == 'spheroid': 561 # Normal value. 562 if sim_index == None: 563 # Fold phi inside 0 and pi. 564 if cdp.diff_tensor.phi >= pi: 565 theta, phi = fold_spherical_angles(cdp.diff_tensor.theta, cdp.diff_tensor.phi) 566 cdp.diff_tensor.set(param='theta', value=theta) 567 cdp.diff_tensor.set(param='phi', value=phi) 568 569 # Simulated theta and phi values. 570 else: 571 # Fold phi_sim inside phi-pi/2 and phi+pi/2. 572 if cdp.diff_tensor.phi_sim[sim_index] >= cdp.diff_tensor.phi + pi/2.0: 573 cdp.diff_tensor.set(param='theta', value=pi - cdp.diff_tensor.theta_sim[sim_index], category='sim', sim_index=sim_index) 574 cdp.diff_tensor.set(param='phi', value=cdp.diff_tensor.phi_sim[sim_index] - pi, category='sim', sim_index=sim_index) 575 elif cdp.diff_tensor.phi_sim[sim_index] <= cdp.diff_tensor.phi - pi/2.0: 576 cdp.diff_tensor.set(param='theta', value=pi - cdp.diff_tensor.theta_sim[sim_index], category='sim', sim_index=sim_index) 577 cdp.diff_tensor.set(param='phi', value=cdp.diff_tensor.phi_sim[sim_index] + pi, category='sim', sim_index=sim_index) 578 579 # Ellipsoid. 580 elif cdp.diff_tensor.type == 'ellipsoid': 581 # Normal value. 582 if sim_index == None: 583 # Fold alpha inside 0 and pi. 584 if cdp.diff_tensor.alpha >= pi: 585 cdp.diff_tensor.set(param='alpha', value=cdp.diff_tensor.alpha - pi) 586 587 # Fold beta inside 0 and pi. 588 if cdp.diff_tensor.beta >= pi: 589 cdp.diff_tensor.set(param='alpha', value=pi - cdp.diff_tensor.alpha) 590 cdp.diff_tensor.set(param='beta', value=cdp.diff_tensor.beta - pi) 591 592 # Fold gamma inside 0 and pi. 593 if cdp.diff_tensor.gamma >= pi: 594 cdp.diff_tensor.set(param='alpha', value=pi - cdp.diff_tensor.alpha) 595 cdp.diff_tensor.set(param='beta', value=pi - cdp.diff_tensor.beta) 596 cdp.diff_tensor.set(param='gamma', value=cdp.diff_tensor.gamma - pi) 597 598 # Simulated theta and phi values. 599 else: 600 # Fold alpha_sim inside alpha-pi/2 and alpha+pi/2. 601 if cdp.diff_tensor.alpha_sim[sim_index] >= cdp.diff_tensor.alpha + pi/2.0: 602 cdp.diff_tensor.set(param='alpha', value=cdp.diff_tensor.alpha_sim[sim_index] - pi, category='sim', sim_index=sim_index) 603 elif cdp.diff_tensor.alpha_sim[sim_index] <= cdp.diff_tensor.alpha - pi/2.0: 604 cdp.diff_tensor.set(param='alpha', value=cdp.diff_tensor.alpha_sim[sim_index] + pi, category='sim', sim_index=sim_index) 605 606 # Fold beta_sim inside beta-pi/2 and beta+pi/2. 607 if cdp.diff_tensor.beta_sim[sim_index] >= cdp.diff_tensor.beta + pi/2.0: 608 cdp.diff_tensor.set(param='alpha', value=pi - cdp.diff_tensor.alpha_sim[sim_index], category='sim', sim_index=sim_index) 609 cdp.diff_tensor.set(param='beta', value=cdp.diff_tensor.beta_sim[sim_index] - pi, category='sim', sim_index=sim_index) 610 elif cdp.diff_tensor.beta_sim[sim_index] <= cdp.diff_tensor.beta - pi/2.0: 611 cdp.diff_tensor.set(param='alpha', value=pi - cdp.diff_tensor.alpha_sim[sim_index], category='sim', sim_index=sim_index) 612 cdp.diff_tensor.set(param='beta', value=cdp.diff_tensor.beta_sim[sim_index] + pi, category='sim', sim_index=sim_index) 613 614 # Fold gamma_sim inside gamma-pi/2 and gamma+pi/2. 615 if cdp.diff_tensor.gamma_sim[sim_index] >= cdp.diff_tensor.gamma + pi/2.0: 616 cdp.diff_tensor.set(param='alpha', value=pi - cdp.diff_tensor.alpha_sim[sim_index], category='sim', sim_index=sim_index) 617 cdp.diff_tensor.set(param='beta', value=pi - cdp.diff_tensor.beta_sim[sim_index], category='sim', sim_index=sim_index) 618 cdp.diff_tensor.set(param='gamma', value=cdp.diff_tensor.gamma_sim[sim_index] - pi, category='sim', sim_index=sim_index) 619 elif cdp.diff_tensor.gamma_sim[sim_index] <= cdp.diff_tensor.gamma - pi/2.0: 620 cdp.diff_tensor.set(param='alpha', value=pi - cdp.diff_tensor.alpha_sim[sim_index], category='sim', sim_index=sim_index) 621 cdp.diff_tensor.set(param='beta', value=pi - cdp.diff_tensor.beta_sim[sim_index], category='sim', sim_index=sim_index) 622 cdp.diff_tensor.set(param='gamma', value=cdp.diff_tensor.gamma_sim[sim_index] + pi, category='sim', sim_index=sim_index)

623 624

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

626 """Function for initialising the diffusion tensor. 627 628 @param params: The diffusion tensor parameters. 629 @type params: float 630 @param time_scale: The correlation time scaling value. 631 @type time_scale: float 632 @param d_scale: The diffusion tensor eigenvalue scaling value. 633 @type d_scale: float 634 @param angle_units: The units for the angle parameters. 635 @type angle_units: str (either 'deg' or 'rad') 636 @param param_types: The type of parameters supplied. For spherical diffusion, the flag 637 values correspond to 0: tm, 1: Diso. For spheroidal diffusion, 0: {tm, 638 Da, theta, phi}, 1: {Diso, Da, theta, phi}, 2: {tm, Dratio, theta, phi}, 639 3: {Dpar, Dper, theta, phi}, 4: {Diso, Dratio, theta, phi}. For 640 ellipsoidal diffusion, 0: {tm, Da, Dr, alpha, beta, gamma}, 1: {Diso, 641 Da, Dr, alpha, beta, gamma}, 2: {Dx, Dy, Dz, alpha, beta, gamma}. 642 @type param_types: int 643 @param spheroid_type: A string which, if supplied together with spheroid parameters, will 644 restrict the tensor to either being 'oblate' or 'prolate'. 645 @type spheroid_type: str 646 @param fixed: A flag specifying whether the diffusion tensor is fixed or can be 647 optimised. 648 @type fixed: bin 649 """ 650 651 # Test if the current data pipe exists. 652 check_pipe() 653 654 # Test if diffusion tensor data already exists. 655 if diff_data_exists(): 656 raise RelaxTensorError('diffusion') 657 658 # Check the validity of the angle_units argument. 659 valid_types = ['deg', 'rad'] 660 if not angle_units in valid_types: 661 raise RelaxError("The diffusion tensor 'angle_units' argument " + repr(angle_units) + " should be either 'deg' or 'rad'.") 662 663 # Add the diff_tensor object to the data pipe. 664 cdp.diff_tensor = DiffTensorData() 665 666 # Set the fixed flag. 667 cdp.diff_tensor.set_fixed(fixed) 668 669 # Spherical diffusion. 670 if isinstance(params, float) or params == None: 671 num_params = 1 672 sphere(params, time_scale, param_types) 673 674 # Spheroidal diffusion. 675 elif (isinstance(params, tuple) or isinstance(params, list)) and len(params) == 4: 676 num_params = 4 677 spheroid(params, time_scale, d_scale, angle_units, param_types, spheroid_type) 678 679 # Ellipsoidal diffusion. 680 elif (isinstance(params, tuple) or isinstance(params, list)) and len(params) == 6: 681 num_params = 6 682 ellipsoid(params, time_scale, d_scale, angle_units, param_types) 683 684 # Unknown. 685 else: 686 raise RelaxError("The diffusion tensor parameters " + repr(params) + " are of an unknown type.") 687 688 # Test the validity of the parameters. 689 test_params(num_params)

690 691

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

693 """The function for creating bounds for the mapping function. 694 695 @param param: The name of the parameter to return the bounds for. 696 @type param: str 697 @keyword spin_id: The spin identification string. This arg is unused. 698 @type spin_id: None or str 699 @return: The bounds for the parameter. 700 @rtype: list of len 2 of floats 701 """ 702 703 # tm. 704 if param == 'tm': 705 return [0, 10.0 * 1e-9] 706 707 # {Diso, Dx, Dy, Dz, Dpar, Dper}. 708 if param == 'Diso' or param == 'Dx' or param == 'Dy' or param == 'Dz' or param == 'Dpar' or param == 'Dper': 709 return [1e6, 1e7] 710 711 # Da. 712 if param == 'Da': 713 return [-3.0/2.0 * 1e7, 3.0 * 1e7] 714 715 # Dr. 716 if param == 'Dr': 717 return [0, 1] 718 719 # Dratio. 720 if param == 'Dratio': 721 return [1.0/3.0, 3.0] 722 723 # theta. 724 if param == 'theta': 725 return [0, pi] 726 727 # phi. 728 if param == 'phi': 729 return [0, 2*pi] 730 731 # alpha. 732 if param == 'alpha': 733 return [0, 2*pi] 734 735 # beta. 736 if param == 'beta': 737 return [0, pi] 738 739 # gamma. 740 if param == 'gamma': 741 return [0, 2*pi]

742 743

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

745 """Function for creating labels, tick locations, and tick values for an OpenDX map.""" 746 747 # Initialise. 748 labels = "{" 749 tick_locations = [] 750 tick_values = [] 751 n = len(params) 752 axis_incs = 5 753 loc_inc = inc / axis_incs 754 755 # Increment over the model parameters. 756 for i in range(n): 757 # Parameter conversion factors. 758 factor = return_conversion_factor(params[swap[i]]) 759 760 # Parameter units. 761 units = return_units(params[swap[i]]) 762 763 # Labels. 764 if units: 765 labels = labels + "\"" + params[swap[i]] + " (" + units + ")\"" 766 else: 767 labels = labels + "\"" + params[swap[i]] + "\"" 768 769 # Tick values. 770 vals = bounds[swap[i], 0] / factor 771 val_inc = (bounds[swap[i], 1] - bounds[swap[i], 0]) / (axis_incs * factor) 772 773 if i < n - 1: 774 labels = labels + " " 775 else: 776 labels = labels + "}" 777 778 # Tick locations. 779 string = "{" 780 val = 0.0 781 for j in range(axis_incs + 1): 782 string = string + " " + repr(val) 783 val = val + loc_inc 784 string = string + " }" 785 tick_locations.append(string) 786 787 # Tick values. 788 string = "{" 789 for j in range(axis_incs + 1): 790 string = string + "\"" + "%.2f" % vals + "\" " 791 vals = vals + val_inc 792 string = string + "}" 793 tick_values.append(string) 794 795 return labels, tick_locations, tick_values

796 797

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

799 """Set the diffusion tensor parameters. 800 801 @keyword tensor: The diffusion tensor object. 802 @type tensor: DiffTensorData instance 803 @keyword value: The list of values to set the parameters to. 804 @type value: list of float 805 @keyword param: The list of parameter names. 806 @type param: list of str 807 """ 808 809 # Set up the diffusion tensor data if it doesn't exist. 810 if not diff_data_exists(): 811 raise RelaxNoTensorError('diffusion') 812 813 # Initialise. 814 geo_params = [] 815 geo_values = [] 816 orient_params = [] 817 orient_values = [] 818 819 # Loop over the parameters. 820 for i in range(len(param)): 821 # Unknown parameter. 822 if not param[i]: 823 raise RelaxUnknownParamError("diffusion tensor", param[i]) 824 825 # Default value. 826 if value[i] == None: 827 value[i] = default_value(param[i]) 828 829 # Geometric parameter. 830 if param[i] in ['tm', 'Diso', 'Da', 'Dratio', 'Dper', 'Dpar', 'Dr', 'Dx', 'Dy', 'Dz']: 831 geo_params.append(param[i]) 832 geo_values.append(value[i]) 833 834 # Orientational parameter. 835 if param[i] in ['theta', 'phi', 'alpha', 'beta', 'gamma']: 836 orient_params.append(param[i]) 837 orient_values.append(value[i]) 838 839 840 # Spherical diffusion. 841 ###################### 842 843 if cdp.diff_tensor.type == 'sphere': 844 # Geometric parameters. 845 ####################### 846 847 # A single geometric parameter. 848 if len(geo_params) == 1: 849 # The single parameter tm. 850 if geo_params[0] == 'tm': 851 cdp.diff_tensor.set(param='tm', value=geo_values[0]) 852 853 # The single parameter Diso. 854 elif geo_params[0] == 'Diso': 855 cdp.diff_tensor.set(param='tm', value=1.0 / (6.0 * geo_values[0])) 856 857 # Cannot set the single parameter. 858 else: 859 raise RelaxError("The geometric diffusion parameter " + repr(geo_params[0]) + " cannot be set.") 860 861 # More than one geometric parameters. 862 elif len(geo_params) > 1: 863 raise RelaxUnknownParamCombError('geometric parameter set', geo_params) 864 865 866 # Orientational parameters. 867 ########################### 868 869 # ??? 870 if len(orient_params): 871 raise RelaxError("For spherical diffusion, the orientation parameters " + repr(orient_params) + " should not exist.") 872 873 874 # Spheroidal diffusion. 875 ####################### 876 877 elif cdp.diff_tensor.type == 'spheroid': 878 # Geometric parameters. 879 ####################### 880 881 # A single geometric parameter. 882 if len(geo_params) == 1: 883 # The single parameter tm. 884 if geo_params[0] == 'tm': 885 cdp.diff_tensor.set(param='tm', value=geo_values[0]) 886 887 # The single parameter Diso. 888 elif geo_params[0] == 'Diso': 889 cdp.diff_tensor.set(param='tm', value=1.0 / (6.0 * geo_values[0])) 890 891 # The single parameter Da. 892 elif geo_params[0] == 'Da': 893 cdp.diff_tensor.set(param='Da', value=geo_values[0]) 894 895 # The single parameter Dratio. 896 elif geo_params[0] == 'Dratio': 897 Dratio = geo_values[0] 898 cdp.diff_tensor.set(param='Da', value=(Dratio - 1.0) / (2.0 * cdp.diff_tensor.tm * (Dratio + 2.0))) 899 900 # Cannot set the single parameter. 901 else: 902 raise RelaxError("The geometric diffusion parameter " + repr(geo_params[0]) + " cannot be set.") 903 904 # Two geometric parameters. 905 elif len(geo_params) == 2: 906 # The geometric parameter set {tm, Da}. 907 if geo_params.count('tm') == 1 and geo_params.count('Da') == 1: 908 # The parameters. 909 tm = geo_values[geo_params.index('tm')] 910 Da = geo_values[geo_params.index('Da')] 911 912 # Set the internal parameter values. 913 cdp.diff_tensor.set(param='tm', value=tm) 914 cdp.diff_tensor.set(param='Da', value=Da) 915 916 # The geometric parameter set {Diso, Da}. 917 elif geo_params.count('Diso') == 1 and geo_params.count('Da') == 1: 918 # The parameters. 919 Diso = geo_values[geo_params.index('Diso')] 920 Da = geo_values[geo_params.index('Da')] 921 922 # Set the internal parameter values. 923 cdp.diff_tensor.set(param='tm', value=1.0 / (6.0 * Diso)) 924 cdp.diff_tensor.set(param='Da', value=Da) 925 926 # The geometric parameter set {tm, Dratio}. 927 elif geo_params.count('tm') == 1 and geo_params.count('Dratio') == 1: 928 # The parameters. 929 tm = geo_values[geo_params.index('tm')] 930 Dratio = geo_values[geo_params.index('Dratio')] 931 932 # Set the internal parameter values. 933 cdp.diff_tensor.set(param='tm', value=tm) 934 cdp.diff_tensor.set(param='Da', value=(Dratio - 1.0) / (2.0 * tm * (Dratio + 2.0))) 935 936 # The geometric parameter set {Dpar, Dper}. 937 elif geo_params.count('Dpar') == 1 and geo_params.count('Dper') == 1: 938 # The parameters. 939 Dpar = geo_values[geo_params.index('Dpar')] 940 Dper = geo_values[geo_params.index('Dper')] 941 942 # Set the internal parameter values. 943 cdp.diff_tensor.set(param='tm', value=1.0 / (2.0 * (Dpar + 2.0*Dper))) 944 cdp.diff_tensor.set(param='Da', value=Dpar - Dper) 945 946 # The geometric parameter set {Diso, Dratio}. 947 elif geo_params.count('Diso') == 1 and geo_params.count('Dratio') == 1: 948 # The parameters. 949 Diso = geo_values[geo_params.index('Diso')] 950 Dratio = geo_values[geo_params.index('Dratio')] 951 952 # Set the internal parameter values. 953 cdp.diff_tensor.set(param='tm', value=1.0 / (6.0 * Diso)) 954 cdp.diff_tensor.set(param='Da', value=3.0 * Diso * (Dratio - 1.0) / (Dratio + 2.0)) 955 956 # Unknown parameter combination. 957 else: 958 raise RelaxUnknownParamCombError('geometric parameter set', geo_params) 959 960 # More than two geometric parameters. 961 elif len(geo_params) > 2: 962 raise RelaxUnknownParamCombError('geometric parameter set', geo_params) 963 964 965 # Orientational parameters. 966 ########################### 967 968 # A single orientational parameter. 969 if len(orient_params) == 1: 970 # The single parameter theta. 971 if orient_params[0] == 'theta': 972 cdp.diff_tensor.set(param='theta', value=orient_values[orient_params.index('theta')]) 973 974 # The single parameter phi. 975 elif orient_params[0] == 'phi': 976 cdp.diff_tensor.set(param='phi', value=orient_values[orient_params.index('phi')]) 977 978 # Disallowed parameter. 979 else: 980 raise RelaxError("For spheroidal diffusion, the orientation parameter " + repr(orient_params) + " cannot be set.") 981 982 # Two orientational parameters. 983 elif len(orient_params) == 2: 984 # The orientational parameter set {theta, phi}. 985 if orient_params.count('theta') == 1 and orient_params.count('phi') == 1: 986 cdp.diff_tensor.set(param='theta', value=orient_values[orient_params.index('theta')]) 987 cdp.diff_tensor.set(param='phi', value=orient_values[orient_params.index('phi')]) 988 989 # Unknown parameter combination. 990 else: 991 raise RelaxUnknownParamCombError('orientational parameter set', orient_params) 992 993 # More than two orientational parameters. 994 elif len(orient_params) > 2: 995 raise RelaxUnknownParamCombError('orientational parameter set', orient_params) 996 997 998 # Ellipsoidal diffusion. 999 ######################## 1000 1001 elif cdp.diff_tensor.type == 'ellipsoid': 1002 # Geometric parameters. 1003 ####################### 1004 1005 # A single geometric parameter. 1006 if len(geo_params) == 1: 1007 # The single parameter tm. 1008 if geo_params[0] == 'tm': 1009 cdp.diff_tensor.set(param='tm', value=geo_values[0]) 1010 1011 # The single parameter Diso. 1012 elif geo_params[0] == 'Diso': 1013 cdp.diff_tensor.set(param='tm', value=1.0 / (6.0 * geo_values[0])) 1014 1015 # The single parameter Da. 1016 elif geo_params[0] == 'Da': 1017 cdp.diff_tensor.set(param='Da', value=geo_values[0]) 1018 1019 # The single parameter Dr. 1020 elif geo_params[0] == 'Dr': 1021 cdp.diff_tensor.set(param='Dr', value=geo_values[0]) 1022 1023 # Cannot set the single parameter. 1024 else: 1025 raise RelaxError("The geometric diffusion parameter " + repr(geo_params[0]) + " cannot be set.") 1026 1027 # Two geometric parameters. 1028 elif len(geo_params) == 2: 1029 # The geometric parameter set {tm, Da}. 1030 if geo_params.count('tm') == 1 and geo_params.count('Da') == 1: 1031 # The parameters. 1032 tm = geo_values[geo_params.index('tm')] 1033 Da = geo_values[geo_params.index('Da')] 1034 1035 # Set the internal parameter values. 1036 cdp.diff_tensor.set(param='tm', value=tm) 1037 cdp.diff_tensor.set(param='Da', value=Da) 1038 1039 # The geometric parameter set {tm, Dr}. 1040 elif geo_params.count('tm') == 1 and geo_params.count('Dr') == 1: 1041 # The parameters. 1042 tm = geo_values[geo_params.index('tm')] 1043 Dr = geo_values[geo_params.index('Dr')] 1044 1045 # Set the internal parameter values. 1046 cdp.diff_tensor.set(param='tm', value=tm) 1047 cdp.diff_tensor.set(param='Dr', value=Dr) 1048 1049 # The geometric parameter set {Diso, Da}. 1050 elif geo_params.count('Diso') == 1 and geo_params.count('Da') == 1: 1051 # The parameters. 1052 Diso = geo_values[geo_params.index('Diso')] 1053 Da = geo_values[geo_params.index('Da')] 1054 1055 # Set the internal parameter values. 1056 cdp.diff_tensor.set(param='tm', value=1.0 / (6.0 * Diso)) 1057 cdp.diff_tensor.set(param='Da', value=Da) 1058 1059 # The geometric parameter set {Diso, Dr}. 1060 elif geo_params.count('Diso') == 1 and geo_params.count('Dr') == 1: 1061 # The parameters. 1062 Diso = geo_values[geo_params.index('Diso')] 1063 Dr = geo_values[geo_params.index('Dr')] 1064 1065 # Set the internal parameter values. 1066 cdp.diff_tensor.set(param='tm', value=1.0 / (6.0 * Diso)) 1067 cdp.diff_tensor.set(param='Dr', value=Dr) 1068 1069 # The geometric parameter set {Da, Dr}. 1070 elif geo_params.count('Da') == 1 and geo_params.count('Dr') == 1: 1071 # The parameters. 1072 Da = geo_values[geo_params.index('Da')] 1073 Dr = geo_values[geo_params.index('Dr')] 1074 1075 # Set the internal parameter values. 1076 cdp.diff_tensor.set(param='Da', value=Da) 1077 cdp.diff_tensor.set(param='Da', value=Dr) 1078 1079 # Unknown parameter combination. 1080 else: 1081 raise RelaxUnknownParamCombError('geometric parameter set', geo_params) 1082 1083 # Three geometric parameters. 1084 elif len(geo_params) == 3: 1085 # The geometric parameter set {tm, Da, Dr}. 1086 if geo_params.count('tm') == 1 and geo_params.count('Da') == 1 and geo_params.count('Dr') == 1: 1087 # The parameters. 1088 tm = geo_values[geo_params.index('tm')] 1089 Da = geo_values[geo_params.index('Da')] 1090 Dr = geo_values[geo_params.index('Dr')] 1091 1092 # Set the internal parameter values. 1093 cdp.diff_tensor.set(param='tm', value=tm) 1094 cdp.diff_tensor.set(param='Da', value=Da) 1095 cdp.diff_tensor.set(param='Dr', value=Dr) 1096 1097 # The geometric parameter set {Diso, Da, Dr}. 1098 elif geo_params.count('Diso') == 1 and geo_params.count('Da') == 1 and geo_params.count('Dr') == 1: 1099 # The parameters. 1100 Diso = geo_values[geo_params.index('Diso')] 1101 Da = geo_values[geo_params.index('Da')] 1102 Dr = geo_values[geo_params.index('Dr')] 1103 1104 # Set the internal parameter values. 1105 cdp.diff_tensor.set(param='tm', value=1.0 / (6.0 * Diso)) 1106 cdp.diff_tensor.set(param='Da', value=Da) 1107 cdp.diff_tensor.set(param='Dr', value=Dr) 1108 1109 # The geometric parameter set {Dx, Dy, Dz}. 1110 elif geo_params.count('Dx') == 1 and geo_params.count('Dy') == 1 and geo_params.count('Dz') == 1: 1111 # The parameters. 1112 Dx = geo_values[geo_params.index('Dx')] 1113 Dy = geo_values[geo_params.index('Dy')] 1114 Dz = geo_values[geo_params.index('Dz')] 1115 1116 # Set the internal tm value. 1117 if Dx + Dy + Dz == 0.0: 1118 cdp.diff_tensor.set(param='tm', value=1e99) 1119 else: 1120 cdp.diff_tensor.set(param='tm', value=0.5 / (Dx + Dy + Dz)) 1121 1122 # Set the internal Da value. 1123 cdp.diff_tensor.set(param='Da', value=Dz - 0.5*(Dx + Dy)) 1124 1125 # Set the internal Dr value. 1126 if cdp.diff_tensor.Da == 0.0: 1127 cdp.diff_tensor.set(param='Dr', value=(Dy - Dx) * 1e99) 1128 else: 1129 cdp.diff_tensor.set(param='Dr', value=(Dy - Dx) / (2.0*cdp.diff_tensor.Da)) 1130 1131 # Unknown parameter combination. 1132 else: 1133 raise RelaxUnknownParamCombError('geometric parameter set', geo_params) 1134 1135 1136 # More than three geometric parameters. 1137 elif len(geo_params) > 3: 1138 raise RelaxUnknownParamCombError('geometric parameter set', geo_params) 1139 1140 1141 # Orientational parameters. 1142 ########################### 1143 1144 # A single orientational parameter. 1145 if len(orient_params) == 1: 1146 # The single parameter alpha. 1147 if orient_params[0] == 'alpha': 1148 cdp.diff_tensor.set(param='alpha', value=orient_values[orient_params.index('alpha')]) 1149 1150 # The single parameter beta. 1151 elif orient_params[0] == 'beta': 1152 cdp.diff_tensor.set(param='beta', value=orient_values[orient_params.index('beta')]) 1153 1154 # The single parameter gamma. 1155 elif orient_params[0] == 'gamma': 1156 cdp.diff_tensor.set(param='gamma', value=orient_values[orient_params.index('gamma')]) 1157 1158 # Disallowed parameter. 1159 else: 1160 raise RelaxError("For ellipsoidal diffusion, the orientation parameter " + repr(orient_params) + " cannot be set.") 1161 1162 # Two orientational parameters. 1163 elif len(orient_params) == 2: 1164 # The orientational parameter set {alpha, beta}. 1165 if orient_params.count('alpha') == 1 and orient_params.count('beta') == 1: 1166 cdp.diff_tensor.set(param='alpha', value=orient_values[orient_params.index('alpha')]) 1167 cdp.diff_tensor.set(param='beta', value=orient_values[orient_params.index('beta')]) 1168 1169 # The orientational parameter set {alpha, gamma}. 1170 if orient_params.count('alpha') == 1 and orient_params.count('gamma') == 1: 1171 cdp.diff_tensor.set(param='alpha', value=orient_values[orient_params.index('alpha')]) 1172 cdp.diff_tensor.set(param='gamma', value=orient_values[orient_params.index('gamma')]) 1173 1174 # The orientational parameter set {beta, gamma}. 1175 if orient_params.count('beta') == 1 and orient_params.count('gamma') == 1: 1176 cdp.diff_tensor.set(param='beta', value=orient_values[orient_params.index('beta')]) 1177 cdp.diff_tensor.set(param='gamma', value=orient_values[orient_params.index('gamma')]) 1178 1179 # Unknown parameter combination. 1180 else: 1181 raise RelaxUnknownParamCombError('orientational parameter set', orient_params) 1182 1183 # Three orientational parameters. 1184 elif len(orient_params) == 3: 1185 # The orientational parameter set {alpha, beta, gamma}. 1186 if orient_params.count('alpha') == 1 and orient_params.count('beta') == 1: 1187 cdp.diff_tensor.set(param='alpha', value=orient_values[orient_params.index('alpha')]) 1188 cdp.diff_tensor.set(param='beta', value=orient_values[orient_params.index('beta')]) 1189 cdp.diff_tensor.set(param='gamma', value=orient_values[orient_params.index('gamma')]) 1190 1191 # Unknown parameter combination. 1192 else: 1193 raise RelaxUnknownParamCombError('orientational parameter set', orient_params) 1194 1195 # More than three orientational parameters. 1196 elif len(orient_params) > 3: 1197 raise RelaxUnknownParamCombError('orientational parameter set', orient_params) 1198 1199 1200 # Fold the angles in. 1201 ##################### 1202 1203 if orient_params: 1204 fold_angles()

1205 1206

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

1208 """Function for setting up a spherical diffusion tensor. 1209 1210 @param params: The diffusion tensor parameter. 1211 @type params: float 1212 @param time_scale: The correlation time scaling value. 1213 @type time_scale: float 1214 @param param_types: The type of parameter supplied. If 0, then the parameter is tm. If 1, then 1215 the parameter is Diso. 1216 @type param_types: int 1217 """ 1218 1219 # The diffusion type. 1220 cdp.diff_tensor.set_type('sphere') 1221 1222 # Set the parameter values. 1223 if params != None: 1224 # tm. 1225 if param_types == 0: 1226 # Scaling. 1227 tm = params * time_scale 1228 1229 # Set the parameters. 1230 set(value=[tm], param=['tm']) 1231 1232 # Diso 1233 elif param_types == 1: 1234 # Scaling. 1235 Diso = params * d_scale 1236 1237 # Set the parameters. 1238 set(value=[Diso], param=['Diso']) 1239 1240 # Unknown parameter combination. 1241 else: 1242 raise RelaxUnknownParamCombError('param_types', param_types)

1243 1244

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

1246 """Function for setting up a spheroidal diffusion tensor. 1247 1248 @param params: The diffusion tensor parameter. 1249 @type params: float 1250 @param time_scale: The correlation time scaling value. 1251 @type time_scale: float 1252 @param d_scale: The diffusion tensor eigenvalue scaling value. 1253 @type d_scale: float 1254 @param angle_units: The units for the angle parameters which can be either 'deg' or 'rad'. 1255 @type angle_units: str 1256 @param param_types: The type of parameters supplied. These correspond to 0: {tm, Da, theta, 1257 phi}, 1: {Diso, Da, theta, phi}, 2: {tm, Dratio, theta, phi}, 3: {Dpar, 1258 Dper, theta, phi}, 4: {Diso, Dratio, theta, phi}. 1259 @type param_types: int 1260 @param spheroid_type: A string which, if supplied together with spheroid parameters, will 1261 restrict the tensor to either being 'oblate' or 'prolate'. 1262 @type spheroid_type: str 1263 """ 1264 1265 # The diffusion type. 1266 cdp.diff_tensor.set_type('spheroid') 1267 1268 # Spheroid diffusion type. 1269 allowed_types = [None, 'oblate', 'prolate'] 1270 if spheroid_type not in allowed_types: 1271 raise RelaxError("The 'spheroid_type' argument " + repr(spheroid_type) + " should be 'oblate', 'prolate', or None.") 1272 cdp.diff_tensor.set(param='spheroid_type', value=spheroid_type) 1273 1274 # Set the parameter values. 1275 if params[0] != None: 1276 # (tm, Da, theta, phi). 1277 if param_types == 0: 1278 # Unpack the tuple. 1279 tm, Da, theta, phi = params 1280 1281 # Scaling. 1282 tm = tm * time_scale 1283 Da = Da * d_scale 1284 1285 # Set the parameters. 1286 set(value=[tm, Da], param=['tm', 'Da']) 1287 1288 # (Diso, Da, theta, phi). 1289 elif param_types == 1: 1290 # Unpack the tuple. 1291 Diso, Da, theta, phi = params 1292 1293 # Scaling. 1294 Diso = Diso * d_scale 1295 Da = Da * d_scale 1296 1297 # Set the parameters. 1298 set(value=[Diso, Da], param=['Diso', 'Da']) 1299 1300 # (tm, Dratio, theta, phi). 1301 elif param_types == 2: 1302 # Unpack the tuple. 1303 tm, Dratio, theta, phi = params 1304 1305 # Scaling. 1306 tm = tm * time_scale 1307 1308 # Set the parameters. 1309 set(value=[tm, Dratio], param=['tm', 'Dratio']) 1310 1311 # (Dpar, Dper, theta, phi). 1312 elif param_types == 3: 1313 # Unpack the tuple. 1314 Dpar, Dper, theta, phi = params 1315 1316 # Scaling. 1317 Dpar = Dpar * d_scale 1318 Dper = Dper * d_scale 1319 1320 # Set the parameters. 1321 set(value=[Dpar, Dper], param=['Dpar', 'Dper']) 1322 1323 # (Diso, Dratio, theta, phi). 1324 elif param_types == 4: 1325 # Unpack the tuple. 1326 Diso, Dratio, theta, phi = params 1327 1328 # Scaling. 1329 Diso = Diso * d_scale 1330 1331 # Set the parameters. 1332 set(value=[Diso, Dratio], param=['Diso', 'Dratio']) 1333 1334 # Unknown parameter combination. 1335 else: 1336 raise RelaxUnknownParamCombError('param_types', param_types) 1337 1338 # Convert the angles to radians. 1339 if angle_units == 'deg': 1340 print("Converting the angles from degrees to radian units.") 1341 theta = (theta / 360.0) * 2.0 * pi 1342 phi = (phi / 360.0) * 2.0 * pi 1343 1344 # Set the orientational parameters. 1345 set(value=[theta, phi], param=['theta', 'phi'])

1346 1347

1348 -def test_params(num_params):

1349 """Function for testing the validity of the input parameters.""" 1350 1351 # No parameters were set. 1352 if not hasattr(cdp.diff_tensor, 'tm'): 1353 return 1354 1355 # An allowable error to account for machine precision, optimisation quality, etc. 1356 error = 1e-4 1357 1358 # tm. 1359 tm = cdp.diff_tensor.tm 1360 if tm <= 0.0 or tm > 1e-6: 1361 raise RelaxError("The tm value of " + repr(tm) + " should be between zero and one microsecond.") 1362 1363 # Spheroid. 1364 if num_params == 4: 1365 # Parameters. 1366 Diso = 1.0 / (6.0 * cdp.diff_tensor.tm) 1367 Da = cdp.diff_tensor.Da 1368 1369 # Da. 1370 if Da < (-1.5*Diso - error*Diso) or Da > (3.0*Diso + error*Diso): 1371 raise RelaxError("The Da value of " + repr(Da) + " should be between -3/2 * Diso and 3Diso.") 1372 1373 # Ellipsoid. 1374 if num_params == 6: 1375 # Parameters. 1376 Diso = 1.0 / (6.0 * cdp.diff_tensor.tm) 1377 Da = cdp.diff_tensor.Da 1378 Dr = cdp.diff_tensor.Dr 1379 1380 # Da. 1381 if Da < (0.0 - error*Diso) or Da > (3.0*Diso + error*Diso): 1382 raise RelaxError("The Da value of " + repr(Da) + " should be between zero and 3Diso.") 1383 1384 # Dr. 1385 if Dr < (0.0 - error) or Dr > (1.0 + error): 1386 raise RelaxError("The Dr value of " + repr(Dr) + " should be between zero and one.")

1387 1388

1389 -def unit_axes():

1390 """Function for calculating the unit axes of the diffusion tensor. 1391 1392 Spheroid 1393 ======== 1394 1395 The unit Dpar vector is:: 1396 1397 | sin(theta) * cos(phi) | 1398 Dpar = | sin(theta) * sin(phi) | 1399 | cos(theta) | 1400 1401 1402 Ellipsoid 1403 ========= 1404 1405 The unit Dx vector is:: 1406 1407 | -sin(alpha) * sin(gamma) + cos(alpha) * cos(beta) * cos(gamma) | 1408 Dx = | -sin(alpha) * cos(gamma) - cos(alpha) * cos(beta) * sin(gamma) | 1409 | cos(alpha) * sin(beta) | 1410 1411 The unit Dy vector is:: 1412 1413 | cos(alpha) * sin(gamma) + sin(alpha) * cos(beta) * cos(gamma) | 1414 Dy = | cos(alpha) * cos(gamma) - sin(alpha) * cos(beta) * sin(gamma) | 1415 | sin(alpha) * sin(beta) | 1416 1417 The unit Dz vector is:: 1418 1419 | -sin(beta) * cos(gamma) | 1420 Dz = | sin(beta) * sin(gamma) | 1421 | cos(beta) | 1422 1423 """ 1424 1425 # Spheroid. 1426 if cdp.diff_tensor.type == 'spheroid': 1427 # Initialise. 1428 Dpar = zeros(3, float64) 1429 1430 # Trig. 1431 sin_theta = sin(cdp.diff_tensor.theta) 1432 cos_theta = cos(cdp.diff_tensor.theta) 1433 sin_phi = sin(cdp.diff_tensor.phi) 1434 cos_phi = cos(cdp.diff_tensor.phi) 1435 1436 # Unit Dpar axis. 1437 Dpar[0] = sin_theta * cos_phi 1438 Dpar[1] = sin_theta * sin_phi 1439 Dpar[2] = cos_theta 1440 1441 # Return the vector. 1442 return Dpar 1443 1444 # Ellipsoid. 1445 if cdp.diff_tensor.type == 'ellipsoid': 1446 # Initialise. 1447 Dx = zeros(3, float64) 1448 Dy = zeros(3, float64) 1449 Dz = zeros(3, float64) 1450 1451 # Trig. 1452 sin_alpha = sin(cdp.diff_tensor.alpha) 1453 cos_alpha = cos(cdp.diff_tensor.alpha) 1454 sin_beta = sin(cdp.diff_tensor.beta) 1455 cos_beta = cos(cdp.diff_tensor.beta) 1456 sin_gamma = sin(cdp.diff_tensor.gamma) 1457 cos_gamma = cos(cdp.diff_tensor.gamma) 1458 1459 # Unit Dx axis. 1460 Dx[0] = -sin_alpha * sin_gamma + cos_alpha * cos_beta * cos_gamma 1461 Dx[1] = -sin_alpha * cos_gamma - cos_alpha * cos_beta * sin_gamma 1462 Dx[2] = cos_alpha * sin_beta 1463 1464 # Unit Dy axis. 1465 Dx[0] = cos_alpha * sin_gamma + sin_alpha * cos_beta * cos_gamma 1466 Dx[1] = cos_alpha * cos_gamma - sin_alpha * cos_beta * sin_gamma 1467 Dx[2] = sin_alpha * sin_beta 1468 1469 # Unit Dz axis. 1470 Dx[0] = -sin_beta * cos_gamma 1471 Dx[1] = sin_beta * sin_gamma 1472 Dx[2] = cos_beta 1473 1474 # Return the vectors. 1475 return Dx, Dy, Dz

1476

Source Code for Module pipe_control.diffusion_tensor