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/rad)", 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