test_suite.shared_data.diffusion_tensor.generate

1 ############################################################################### 2 # # 3 # Copyright (C) 2010 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 """relax script for creating a PDB file and relaxation data. 23 24 The PDB file consists of uniformly distributed bond vectors. The relaxation data is that of a NH bond vector with an ellipsoidal diffusion tensor and no internal motion. 25 """ 26 27 # Python module imports. 28 from math import pi, sqrt 29 from numpy import array, cross, dot, eye, float64, transpose, zeros 30 from numpy.linalg import eig, inv, norm 31 32 # relax module imports. 33 from maths_fns.rotation_matrix import axis_angle_to_R, euler_to_R_zyz, R_to_euler_zyz 34 from generic_fns.structure.geometric import angles_uniform, vect_dist_spherical_angles 35 from generic_fns.structure.internal import Internal 36 from relax_io import open_write_file 37 38

39 -def ri_data(Dx=None, Dy=None, Dz=None, R=eye(3), vectors=None, frq_label=None, wH=None, csa=None):

40 """Calculate the relaxation data for the given vectors.""" 41 42 # Diff parameters. 43 Diso = (Dx + Dy + Dz) / 3.0 44 L2 = (Dx*Dy + Dx*Dz + Dy*Dz) / 3.0 45 fact = sqrt(Diso**2 - L2) 46 if fact == 0.0: 47 mux = muy = muz = 0.0 48 else: 49 mux = (Dx - Diso) / fact 50 muy = (Dy - Diso) / fact 51 muz = (Dz - Diso) / fact 52 53 # The five time constants. 54 tau = zeros(5, float64) 55 tau[0] = 6 * (Diso - sqrt(Diso**2 - L2)) 56 tau[1] = 4*Dx + Dy + Dz 57 tau[2] = Dx + 4*Dy + Dz 58 tau[3] = Dx + Dy + 4*Dz 59 tau[4] = 6 * (Diso + sqrt(Diso**2 - L2)) 60 tau = 1.0 / tau 61 print("\nTime constants, tau: %s" % tau) 62 63 # The dipolar constant. 64 h = 6.62606876e-34 # Planck constant. 65 h_bar = h / ( 2.0*pi ) # Dirac constant. 66 mu0 = 4.0 * pi * 1e-7 # Permeability of free space. 67 r = 1.02e-10 # NH bond length. 68 gn = -2.7126e7 # 15N gyromagnetic ratio. 69 gh = 26.7522212e7 # 1H gyromagnetic ratio. 70 dip_const = 0.25 * (mu0/(4.0*pi))**2 * (gn * gh * h_bar)**2 / r**6 # The dipolar constant. 71 print("Dipolar constant: %s" % dip_const) 72 73 # The five frequencies. 74 wN = wH * gn/gh 75 w = zeros(5, float64) 76 w[0] = 0 77 w[1] = wN 78 w[2] = wH - wN 79 w[3] = wH 80 w[4] = wH + wN 81 print("Frequencies, w: %s" % w) 82 83 # CSA constant. 84 csa_const = (wN * csa)**2 / 3.0 85 print("CSA constant: %s" % csa_const) 86 87 # The files. 88 R1_file = open('R1.%s.out' % frq_label, 'w') 89 R2_file = open('R2.%s.out' % frq_label, 'w') 90 NOE_file = open('NOE.%s.out' % frq_label, 'w') 91 92 # Loop over the vectors. 93 c = zeros(5, float64) 94 for i in range(len(vectors)): 95 # Normalise. 96 vector = vectors[i] / norm(vectors[i]) 97 98 # Rotate into the diffusion frame. 99 vector = dot(R, vector) 100 101 # Print out. 102 print("\ni: %s" % i) 103 print("vector: %s" % vector) 104 105 # Direction cosines. 106 delta_x, delta_y, delta_z = vector 107 108 # The d and e factors. 109 d = 3.0 * (delta_x**4 + delta_y**4 + delta_z**4) - 1.0 110 e = mux * (delta_x**4 + 2.0*delta_y**2 * delta_z**2) 111 e = e + muy * (delta_y**4 + 2.0*delta_x**2 * delta_z**2) 112 e = e + muz * (delta_z**4 + 2.0*delta_x**2 * delta_y**2) 113 114 # The weights. 115 c[0] = (d + e) / 4.0 116 c[1] = 3 * delta_y**2 * delta_z**2 117 c[2] = 3 * delta_x**2 * delta_z**2 118 c[3] = 3 * delta_x**2 * delta_y**2 119 c[4] = (d - e) / 4.0 120 print("Weights, c: %s" % c) 121 122 # The spectral density function. 123 Jw = zeros(5, float64) 124 for frq_index in range(5): 125 for k in range(5): 126 Jw[frq_index] = Jw[frq_index] + 0.4 * c[k] * tau[k] / (1.0 + (w[frq_index]*tau[k])**2) 127 print("Jw: %s" % Jw) 128 129 # The relaxation data. 130 R1 = dip_const * (Jw[2] + 3.0*Jw[1] + 6.0*Jw[4]) + csa_const * Jw[1] 131 R2 = dip_const/2.0 * (4.0*Jw[0] + Jw[2] + 3.0*Jw[1] + 6.0*Jw[3] + 6.0*Jw[4]) + csa_const/6.0 * (4.0*Jw[0] + 3.0*Jw[1]) 132 sigma_noe = dip_const * (6.0*Jw[4] - Jw[2]) 133 NOE = 1.0 + gh/gn * sigma_noe / R1 134 print("R1: %s" % R1) 135 print("R2: %s" % R2) 136 print("NOE: %s" % NOE) 137 138 # Write the data. 139 R1_file.write("%s %s %s\n" % (i+1, R1, 0.05*R1)) 140 R2_file.write("%s %s %s\n" % (i+1, R2, 0.05*R2)) 141 NOE_file.write("%s %s %s\n" % (i+1, NOE, 0.04))

142 143

144 -def tensor_setup(Dx=None, Dy=None, Dz=None, alpha=None, beta=None, gamma=None):

145 """Set up the diffusion tensor according to the correct Euler angle convention.""" 146 147 # Print out. 148 print("\n\n") 149 print("# Angles to diff tensor.") 150 print("########################") 151 152 # Init. 153 ROT = False 154 155 # The rotation matrix (in the rotating axis system). 156 R = zeros((3, 3), float64) 157 R_rev = zeros((3, 3), float64) 158 euler_to_R_zyz(gamma, beta, alpha, R) 159 R_rev = transpose(R) 160 print("\nEuler angels: [%s, %s, %s]" % (alpha, beta, gamma)) 161 print("R:\n%s" % R) 162 print("R_rev:\n%s" % R_rev) 163 print("X x Y: %s" % cross(R[:, 0], R[:, 1])) 164 165 # Axis rotations. 166 if ROT: 167 R_x180 = zeros((3, 3), float64) 168 R_y180 = zeros((3, 3), float64) 169 R_z180 = zeros((3, 3), float64) 170 axis_angle_to_R(R_rev[:, 0], pi, R_x180) 171 axis_angle_to_R(R_rev[:, 1], pi, R_y180) 172 axis_angle_to_R(R_rev[:, 2], pi, R_z180) 173 print("\nR (x 180):\n%s" % R_x180) 174 print("\nR (y 180):\n%s" % R_y180) 175 print("\nR (z 180):\n%s" % R_z180) 176 177 # A test vector. 178 mu = array([1, 2, -3], float64) 179 mu = mu / norm(mu) 180 181 # Tensor in eigenframe. 182 D_prime = zeros((3, 3), float64) 183 D_prime[0, 0] = Dx 184 D_prime[1, 1] = Dy 185 D_prime[2, 2] = Dz 186 print("\nD':\n%s" % D_prime) 187 188 # Rotate tensor from the eigenframe to the ref frame. 189 D = dot(R_rev, dot(D_prime, transpose(R_rev))) 190 print("\nD:\n%s" % D) 191 print("\n\n") 192 193 # Return the forward and reverse rotation, and the diffusion tensors. 194 return R, R_rev, D_prime, D

195 196

197 -def pdb(r=1.02, file_name='uniform.pdb', inc=None):

198 """Create the bond vector distribution and save the PDB file.""" 199 200 # Create the structural object. 201 structure = Internal() 202 203 # Add a molecule. 204 structure.add_molecule(name='dist') 205 206 # Alias the single molecule from the single model. 207 mol = structure.structural_data[0].mol[0] 208 209 # Get the polar and azimuthal angles for the distribution. 210 phi, theta = angles_uniform(inc) 211 212 # Get the uniform vector distribution. 213 vectors = vect_dist_spherical_angles(inc=inc, distribution='uniform') 214 215 # Loop over the radial array of vectors (change in longitude). 216 atom_num = 1 217 new_vectors = [] 218 for i in range(len(theta)): 219 # Loop over the vectors of the radial array (change in latitude). 220 for j in range(len(phi)): 221 # The index. 222 index = i + j*len(theta) 223 224 # Scale the vector. 225 vector = vectors[index] * r 226 227 # Store the rearranged vector (truncated as in the PDB). 228 trunc_vect = zeros(3, float64) 229 for k in range(3): 230 trunc_vect[k] = float("%.3f" % vector[k]) 231 new_vectors.append(trunc_vect) 232 233 # Residue number. 234 res = (atom_num + 1) / 2 235 236 # Add the vector as a N-H atom pair. 237 mol.atom_add(pdb_record='HETATM', atom_num=atom_num, atom_name='N', res_name='NH', res_num=res, pos=zeros(3), element='N') 238 mol.atom_add(pdb_record='HETATM', atom_num=atom_num+1, atom_name='H', res_name='NH', res_num=res, pos=vector, element='H') 239 240 # Connect. 241 mol.atom_connect(atom_num-1, atom_num) 242 243 # Move 2 atoms forwards. 244 atom_num += 2 245 246 # The PDB file. 247 file = open_write_file(file_name, force=True) 248 structure.write_pdb(file) 249 file.close() 250 251 # Return the vectors in the diffusion frame. 252 return new_vectors

253

Source Code for Module test_suite.shared_data.diffusion_tensor.generate_data