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 41

42 -def bmrb_read(star):

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

139 140

141 -def bmrb_write(star):

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

265 266 267

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

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

305 306

307 -def data_names():

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

314 315

316 -def delete():

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

328 329

330 -def diff_data_exists(pipe=None):

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

351 352

353 -def display():

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

370 371

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

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

468 469

470 -def fold_angles(sim_index=None):

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

620 621

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

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

687 688

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

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

739 740

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

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

793 794

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

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

1202 1203

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

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

1238 1239

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

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

1339 1340

1341 -def test_params(num_params):

1342 """Function for testing the validity of the input parameters.""" 1343 1344 # An allowable error to account for machine precision, optimisation quality, etc. 1345 error = 1e-4 1346 1347 # tm. 1348 tm = cdp.diff_tensor.tm 1349 if tm <= 0.0 or tm > 1e-6: 1350 raise RelaxError("The tm value of " + repr(tm) + " should be between zero and one microsecond.") 1351 1352 # Spheroid. 1353 if num_params == 4: 1354 # Parameters. 1355 Diso = 1.0 / (6.0 * cdp.diff_tensor.tm) 1356 Da = cdp.diff_tensor.Da 1357 1358 # Da. 1359 if Da < (-1.5*Diso - error*Diso) or Da > (3.0*Diso + error*Diso): 1360 raise RelaxError("The Da value of " + repr(Da) + " should be between -3/2 * Diso and 3Diso.") 1361 1362 # Ellipsoid. 1363 if num_params == 6: 1364 # Parameters. 1365 Diso = 1.0 / (6.0 * cdp.diff_tensor.tm) 1366 Da = cdp.diff_tensor.Da 1367 Dr = cdp.diff_tensor.Dr 1368 1369 # Da. 1370 if Da < (0.0 - error*Diso) or Da > (3.0*Diso + error*Diso): 1371 raise RelaxError("The Da value of " + repr(Da) + " should be between zero and 3Diso.") 1372 1373 # Dr. 1374 if Dr < (0.0 - error) or Dr > (1.0 + error): 1375 raise RelaxError("The Dr value of " + repr(Dr) + " should be between zero and one.")

1376 1377

1378 -def unit_axes():

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

1465

Source Code for Module pipe_control.diffusion_tensor