Package lib :: Package diffusion :: Module main
[hide private]
[frames] | no frames]

Source Code for Module lib.diffusion.main

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2003-2008,2010-2011,2013 Edward d'Auvergne                    # 
  4  #                                                                             # 
  5  # This file is part of the program relax (http://www.nmr-relax.com).          # 
  6  #                                                                             # 
  7  # This program is free software: you can redistribute it and/or modify        # 
  8  # it under the terms of the GNU General Public License as published by        # 
  9  # the Free Software Foundation, either version 3 of the License, or           # 
 10  # (at your option) any later version.                                         # 
 11  #                                                                             # 
 12  # This program is distributed in the hope that it will be useful,             # 
 13  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
 14  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
 15  # GNU General Public License for more details.                                # 
 16  #                                                                             # 
 17  # You should have received a copy of the GNU General Public License           # 
 18  # along with this program.  If not, see <http://www.gnu.org/licenses/>.       # 
 19  #                                                                             # 
 20  ############################################################################### 
 21   
 22  # Module docstring. 
 23  """Module for the support of diffusion tensors.""" 
 24   
 25  # Python module imports. 
 26  from math import pi 
 27  from numpy import cross, float64, transpose, zeros 
 28  from numpy.linalg import norm, svd 
 29   
 30  # relax module imports. 
 31  from lib.geometry.rotations import R_to_euler_zyz 
 32  from lib.text.table import format_table 
 33   
 34   
 35   
36 -def return_eigenvalues():
37 """Function for returning Dx, Dy, and Dz.""" 38 39 # Reassign the data. 40 data = cdp.diff_tensor 41 42 # Diso. 43 Diso = 1.0 / (6.0 * data.tm) 44 45 # Dx. 46 Dx = Diso - 1.0/3.0 * data.Da * (1.0 + 3.0 * data.Dr) 47 48 # Dy. 49 Dy = Diso - 1.0/3.0 * data.Da * (1.0 - 3.0 * data.Dr) 50 51 # Dz. 52 Dz = Diso + 2.0/3.0 * data.Da 53 54 # Return the eigenvalues. 55 return Dx, Dy, Dz
56 57
58 -def tensor_eigen_system(tensor):
59 """Determine the eigenvalues and vectors for the tensor, sorting the entries. 60 61 @return: The eigenvalues, rotation matrix, and the Euler angles in zyz notation. 62 @rtype: 3D rank-1 array, 3D rank-2 array, float, float, float 63 """ 64 65 # Eigenvalues. 66 R, Di, A = svd(tensor) 67 D_diag = zeros((3, 3), float64) 68 for i in range(3): 69 D_diag[i, i] = Di[i] 70 71 # Reordering structure. 72 reorder_data = [] 73 for i in range(3): 74 reorder_data.append([Di[i], i]) 75 reorder_data.sort() 76 77 # The indices. 78 reorder = zeros(3, int) 79 Di_sort = zeros(3, float) 80 for i in range(3): 81 Di_sort[i], reorder[i] = reorder_data[i] 82 83 # Reorder columns. 84 R_new = zeros((3, 3), float64) 85 for i in range(3): 86 R_new[:, i] = R[:, reorder[i]] 87 88 # Switch from the left handed to right handed universes (if needed). 89 if norm(cross(R_new[:, 0], R_new[:, 1]) - R_new[:, 2]) > 1e-7: 90 R_new[:, 2] = -R_new[:, 2] 91 92 # Reverse the rotation. 93 R_new = transpose(R_new) 94 95 # Euler angles (reverse rotation in the rotated axis system). 96 gamma, beta, alpha = R_to_euler_zyz(R_new) 97 98 # Collapse the pi axis rotation symmetries. 99 if alpha >= pi: 100 alpha = alpha - pi 101 if gamma >= pi: 102 alpha = pi - alpha 103 beta = pi - beta 104 gamma = gamma - pi 105 if beta >= pi: 106 alpha = pi - alpha 107 beta = beta - pi 108 109 # Return the values. 110 return Di_sort, R_new, alpha, beta, gamma
111 112
113 -def tensor_info_table(type=None, tm=None, Diso=None, Da=None, Dpar=None, Dper=None, Dratio=None, Dr=None, Dx=None, Dy=None, Dz=None, theta=None, phi=None, alpha=None, beta=None, gamma=None, fixed=None):
114 """Print out details of the diffusion tensor. 115 116 @keyword type: The diffusion tensor type - one of 'sphere', 'spheroid', or 'ellipsoid'. 117 @type type: str 118 @keyword tm: The isotropic correlation time in seconds. 119 @type tm: float 120 @keyword Diso: The isotropic diffusion rate. 121 @type Diso: float 122 @keyword Da: The anisotropic component of the tensor. 123 @type Da: float or None 124 @keyword Dpar: The parallel component of the spheroidal diffusion tensor. 125 @type Dpar: float or None 126 @keyword Dper: The perpendicular component of the spheroidal diffusion tensor. 127 @type Dper: float or None 128 @keyword Dratio: The ratio of Dpar and Dper. 129 @type Dratio: float or None 130 @keyword Dr: The rhombic component of the diffusion tensor. 131 @type Dr: float or None 132 @keyword Dx: The x component of the ellipsoid. 133 @type Dx: float or None 134 @keyword Dy: The y component of the ellipsoid. 135 @type Dy: float or None 136 @keyword Dz: The z component of the ellipsoid. 137 @type Dz: float or None 138 @keyword theta: The azimuthal angle in radians. 139 @type theta: float or None 140 @keyword phi: The polar angle in radians. 141 @type phi: float or None 142 @keyword alpha: The Euler angle alpha in radians using the z-y-z convention. 143 @type alpha: float or None 144 @keyword beta: The Euler angle beta in radians using the z-y-z convention. 145 @type beta: float or None 146 @keyword gamma: The Euler angle gamma in radians using the z-y-z convention. 147 @type gamma: float or None 148 """ 149 150 # Build the data for a table. 151 contents = [["Diffusion type", type]] 152 contents.append(["tm (s)", tm]) 153 contents.append(["Diso (rad/s)", Diso]) 154 if Da != None: 155 contents.append(["Da (rad/s)", Da]) 156 if Dpar != None: 157 contents.append(["Dpar (rad/s)", Dpar]) 158 if Dper != None: 159 contents.append(["Dper (rad/s)", Dper]) 160 if Dratio != None: 161 contents.append(["Dratio", Dratio]) 162 if Dr != None: 163 contents.append(["Dr", Dr]) 164 if Dx != None: 165 contents.append(["Dx (rad/s)", Dx]) 166 if Dy != None: 167 contents.append(["Dy (rad/s)", Dy]) 168 if Dz != None: 169 contents.append(["Dz (rad/s)", Dz]) 170 if theta != None: 171 contents.append(["theta (rad)", theta]) 172 if phi != None: 173 contents.append(["phi (rad)", phi]) 174 if alpha != None: 175 contents.append(["alpha (rad)", alpha]) 176 if beta != None: 177 contents.append(["beta (rad)", beta]) 178 if gamma != None: 179 contents.append(["gamma (rad)", gamma]) 180 if fixed != None: 181 contents.append(["Fixed flag", fixed]) 182 183 # Print out the table. 184 print(format_table(contents=contents))
185