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. # 6 # # 7 # relax 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 2 of the License, or # 10 # (at your option) any later version. # 11 # # 12 # relax 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 relax; if not, write to the Free Software # 19 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA # 20 # # 21 ############################################################################### 22 23 """relax script for creating a PDB file and relaxation data. 24 25 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. 26 """ 27 28 # Python module imports. 29 from math import pi, sqrt 30 from numpy import array, cross, dot, eye, float64, transpose, zeros 31 from numpy.linalg import eig, inv, norm 32 33 # relax module imports. 34 from maths_fns.rotation_matrix import axis_angle_to_R, euler_to_R_zyz, R_to_euler_zyz 35 from generic_fns.structure.geometric import angles_uniform, vect_dist_spherical_angles 36 from generic_fns.structure.internal import Internal 37 from relax_io import open_write_file 38 39

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

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

143 144

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

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

196 197

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

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

254

Source Code for Module test_suite.shared_data.diffusion_tensor.generate_data