Package generic_fns :: Module dasha
[hide private]
[frames] | no frames]

Source Code for Module generic_fns.dasha

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2005-2006 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  from math import pi 
 24  from os import F_OK, access, chdir, getcwd, system 
 25  from re import match, search 
 26  from string import lower, split 
 27  import sys 
 28   
 29   
30 -class Dasha:
31 - def __init__(self, relax):
32 """Class used to create and process input and output for the program Modelfree 4.""" 33 34 self.relax = relax
35 36
37 - def create(self, run=None, algor='LM', dir=None, force=0):
38 """Function for creating the Dasha script file 'dir/dasha_script'.""" 39 40 # Arguments. 41 self.run = run 42 self.algor = algor 43 self.dir = dir 44 self.force = force 45 46 # Test if the run exists. 47 if not self.run in self.relax.data.run_names: 48 raise RelaxNoRunError, self.run 49 50 # Test if sequence data is loaded. 51 if not self.relax.data.res.has_key(self.run): 52 raise RelaxNoSequenceError, self.run 53 54 # Determine the parameter set. 55 self.param_set = self.relax.specific.model_free.determine_param_set_type(self.run) 56 57 # Test if diffusion tensor data for the run exists. 58 if self.param_set != 'local_tm' and not self.relax.data.diff.has_key(self.run): 59 raise RelaxNoTensorError, self.run 60 61 # Test if the PDB file has been loaded (for the spheroid and ellipsoid). 62 if self.param_set != 'local_tm' and self.relax.data.diff[self.run].type != 'sphere' and not self.relax.data.pdb.has_key(self.run): 63 raise RelaxNoPdbError, self.run 64 65 # Test if the nucleus type has been set. 66 if not hasattr(self.relax.data, 'gx'): 67 raise RelaxNucleusError 68 69 # Test the optimisation algorithm. 70 if algor not in ['LM', 'NR']: 71 raise RelaxError, "The Dasha optimisation algorithm " + `algor` + " is unknown, it should either be 'LM' or 'NR'." 72 73 # Directory creation. 74 if self.dir == None: 75 self.dir = self.run 76 self.relax.IO.mkdir(self.dir, print_flag=0) 77 78 # Number of field strengths and values. 79 self.num_frq = 0 80 self.frq = [] 81 for i in xrange(len(self.relax.data.res[self.run])): 82 if hasattr(self.relax.data.res[self.run][i], 'num_frq'): 83 if self.relax.data.res[self.run][i].num_frq > self.num_frq: 84 # Number of field strengths. 85 self.num_frq = self.relax.data.res[self.run][i].num_frq 86 87 # Field strength values. 88 for frq in self.relax.data.res[self.run][i].frq: 89 if frq not in self.frq: 90 self.frq.append(frq) 91 92 # Calculate the angle alpha of the XH vector in the spheroid diffusion frame. 93 if self.relax.data.diff[self.run].type == 'spheroid': 94 self.relax.generic.angles.spheroid_frame(self.run) 95 96 # Calculate the angles theta and phi of the XH vector in the ellipsoid diffusion frame. 97 elif self.relax.data.diff[self.run].type == 'ellipsoid': 98 self.relax.generic.angles.ellipsoid_frame(self.run) 99 100 # The 'dasha_script' file. 101 script = self.relax.IO.open_write_file(file_name='dasha_script', dir=self.dir, force=self.force) 102 self.create_script(script) 103 script.close()
104 105
106 - def create_script(self, file):
107 """Create the Dasha script file.""" 108 109 # Delete all data. 110 file.write('# Delete all data.\n') 111 file.write('del 1 10000\n') 112 113 # Nucleus type. 114 file.write('\n# Nucleus type.\n') 115 nucleus = self.relax.generic.nuclei.find_nucleus() 116 if nucleus == 'N': 117 nucleus = 'N15' 118 elif nucleus == 'C': 119 nucleus = 'C13' 120 else: 121 raise RelaxError, 'Cannot handle the nucleus type ' + `nucleus` + ' within Dasha.' 122 file.write('set nucl ' + nucleus + '\n') 123 124 # Number of frequencies. 125 file.write('\n# Number of frequencies.\n') 126 file.write('set n_freq ' + `self.relax.data.num_frq[self.run]` + '\n') 127 128 # Frequency values. 129 file.write('\n# Frequency values.\n') 130 for i in xrange(self.relax.data.num_frq[self.run]): 131 file.write('set H1_freq ' + `self.relax.data.frq[self.run][i] / 1e6` + ' ' + `i+1` + '\n') 132 133 # Set the diffusion tensor. 134 file.write('\n# Set the diffusion tensor.\n') 135 if self.param_set != 'local_tm': 136 # Sphere. 137 if self.relax.data.diff[self.run].type == 'sphere': 138 file.write('set tr ' + `self.relax.data.diff[self.run].tm / 1e-9` + '\n') 139 140 # Spheroid. 141 elif self.relax.data.diff[self.run].type == 'spheroid': 142 file.write('set tr ' + `self.relax.data.diff[self.run].tm / 1e-9` + '\n') 143 144 # Ellipsoid. 145 elif self.relax.data.diff[self.run].type == 'ellipsoid': 146 # Get the eigenvales. 147 Dx, Dy, Dz = self.relax.generic.diffusion_tensor.return_eigenvalues(self.run) 148 149 # Geometric parameters. 150 file.write('set tr ' + `self.relax.data.diff[self.run].tm / 1e-9` + '\n') 151 file.write('set D1/D3 ' + `Dx / Dz` + '\n') 152 file.write('set D2/D3 ' + `Dy / Dz` + '\n') 153 154 # Orientational parameters. 155 file.write('set alfa ' + `self.relax.data.diff[self.run].alpha / (2.0 * pi) * 360.0` + '\n') 156 file.write('set betta ' + `self.relax.data.diff[self.run].beta / (2.0 * pi) * 360.0` + '\n') 157 file.write('set gamma ' + `self.relax.data.diff[self.run].gamma / (2.0 * pi) * 360.0` + '\n') 158 159 # Reading the relaxation data. 160 file.write('\n# Reading the relaxation data.\n') 161 file.write('echo Reading the relaxation data.\n') 162 noe_index = 1 163 r1_index = 1 164 r2_index = 1 165 for i in xrange(self.relax.data.num_ri[self.run]): 166 # NOE. 167 if self.relax.data.ri_labels[self.run][i] == 'NOE': 168 # Data set number. 169 number = noe_index 170 171 # Data type. 172 data_type = 'noe' 173 174 # Increment the data set index. 175 noe_index = noe_index + 1 176 177 # R1. 178 elif self.relax.data.ri_labels[self.run][i] == 'R1': 179 # Data set number. 180 number = r1_index 181 182 # Data type. 183 data_type = '1/T1' 184 185 # Increment the data set index. 186 r1_index = r1_index + 1 187 188 # R2. 189 elif self.relax.data.ri_labels[self.run][i] == 'R2': 190 # Data set number. 191 number = r2_index 192 193 # Data type. 194 data_type = '1/T2' 195 196 # Increment the data set index. 197 r2_index = r2_index + 1 198 199 # Set the data type. 200 if number == 1: 201 file.write('\nread < ' + data_type + '\n') 202 else: 203 file.write('\nread < ' + data_type + ' ' + `number` + '\n') 204 205 # The relaxation data. 206 for j in xrange(len(self.relax.data.res[self.run])): 207 # Reassign the data. 208 data = self.relax.data.res[self.run][j] 209 210 # Skip unselected residues. 211 if not data.select: 212 continue 213 214 # Data and errors. 215 file.write(`data.num` + ' ' + `data.relax_data[i]` + ' ' + `data.relax_error[i]` + '\n') 216 217 # Terminate the reading. 218 file.write('exit\n') 219 220 # Individual residue optimisation. 221 if self.param_set == 'mf': 222 # Loop over the residues. 223 for i in xrange(len(self.relax.data.res[self.run])): 224 # Reassign the data. 225 data = self.relax.data.res[self.run][i] 226 227 # Skip unselected residues. 228 if not data.select: 229 continue 230 231 # Comment. 232 file.write('\n\n\n# Residue ' + `data.num` + '\n\n') 233 234 # Echo. 235 file.write('echo Optimisation of residue ' + `data.num` + '\n') 236 237 # Select the residue. 238 file.write('\n# Select the residue.\n') 239 file.write('set cres ' + `data.num` + '\n') 240 241 # The angle alpha of the XH vector in the spheroid diffusion frame. 242 if self.relax.data.diff[self.run].type == 'spheroid': 243 file.write('set teta ' + `data.alpha` + '\n') 244 245 # The angles theta and phi of the XH vector in the ellipsoid diffusion frame. 246 elif self.relax.data.diff[self.run].type == 'ellipsoid': 247 file.write('\n# Setting the spherical angles of the XH vector in the ellipsoid diffusion frame.\n') 248 file.write('set teta ' + `data.theta` + '\n') 249 file.write('set fi ' + `data.phi` + '\n') 250 251 # The 'jmode'. 252 if 'ts' in data.params: 253 jmode = 3 254 elif 'te' in data.params: 255 jmode = 2 256 elif 'S2' in data.params: 257 jmode = 1 258 259 # Chemical exchange. 260 if 'Rex' in data.params: 261 exch = 1 262 else: 263 exch = 0 264 265 # Anisotropic diffusion. 266 if self.relax.data.diff[self.run].type == 'sphere': 267 anis = 0 268 else: 269 anis = 1 270 271 # Axial symmetry. 272 if self.relax.data.diff[self.run].type == 'spheroid': 273 sym = 1 274 else: 275 sym = 0 276 277 # Set the jmode. 278 file.write('\n# Set the jmode.\n') 279 file.write('set def jmode ' + `jmode`) 280 if exch: 281 file.write(' exch') 282 if anis: 283 file.write(' anis') 284 if sym: 285 file.write(' sym') 286 file.write('\n') 287 288 # Parameter default values. 289 file.write('\n# Parameter default values.\n') 290 file.write('reset jmode ' + `data.num` + '\n') 291 292 # Bond length. 293 file.write('\n# Bond length.\n') 294 file.write('set r_hx ' + `data.r / 1e-10` + '\n') 295 296 # CSA value. 297 file.write('\n# CSA value.\n') 298 file.write('set csa ' + `data.csa / 1e-6` + '\n') 299 300 # Fix the tf parameter if it isn't in the model. 301 if not 'tf' in data.params and jmode == 3: 302 file.write('\n# Fix the tf parameter.\n') 303 file.write('fix tf 0\n') 304 305 # Optimisation of all residues. 306 file.write('\n\n\n# Optimisation of all residues.\n') 307 if self.algor == 'LM': 308 file.write('lmin ' + `self.relax.data.res[self.run][0].num` + ' ' + `self.relax.data.res[self.run][-1].num`) 309 elif self.algor == 'NR': 310 file.write('min ' + `self.relax.data.res[self.run][0].num` + ' ' + `self.relax.data.res[self.run][-1].num`) 311 312 # Show the results. 313 file.write('\n# Show the results.\n') 314 file.write('echo\n') 315 file.write('show all\n') 316 317 # Write the results. 318 file.write('\n# Write the results.\n') 319 file.write('write S2.out S\n') 320 file.write('write S2f.out Sf\n') 321 file.write('write S2s.out Ss\n') 322 file.write('write te.out te\n') 323 file.write('write tf.out tf\n') 324 file.write('write ts.out ts\n') 325 file.write('write Rex.out rex\n') 326 file.write('write chi2.out F\n') 327 328 else: 329 raise RelaxError, 'Optimisation of the parameter set ' + `self.param_set` + ' currently not supported.'
330 331
332 - def execute(self, run, dir, force, binary):
333 """Function for executing Dasha.""" 334 335 # Arguments. 336 self.run = run 337 self.dir = dir 338 self.force = force 339 self.binary = binary 340 341 # Test the binary file string corresponds to a valid executable. 342 self.relax.IO.test_binary(self.binary) 343 344 # The current directory. 345 orig_dir = getcwd() 346 347 # The directory. 348 if self.dir == None: 349 self.dir = self.run 350 if not access(self.dir, F_OK): 351 raise RelaxDirError, ('Dasha', self.dir) 352 353 # Change to this directory. 354 chdir(self.dir) 355 356 # Catch failures and return to the correct directory. 357 try: 358 # Test if the 'dasha_script' script file exists. 359 if not access('dasha_script', F_OK): 360 raise RelaxFileError, ('dasha script', 'dasha_script') 361 362 # Execute Dasha. 363 system(binary + ' < dasha_script | tee dasha_results') 364 365 # Failure. 366 except: 367 # Change back to the original directory. 368 chdir(orig_dir) 369 370 # Reraise the error. 371 raise 372 373 # Change back to the original directory. 374 chdir(orig_dir) 375 376 # Print some blank lines (aesthetics) 377 sys.stdout.write('\n\n')
378 379
380 - def extract(self, run, dir):
381 """Function for extracting the Dasha results out of the 'dasha_results' file.""" 382 383 # Arguments. 384 self.run = run 385 386 # Test if sequence data is loaded. 387 if not self.relax.data.res.has_key(self.run): 388 raise RelaxNoSequenceError, self.run 389 390 # The directory. 391 if dir == None: 392 dir = run 393 if not access(dir, F_OK): 394 raise RelaxDirError, ('Dasha', dir) 395 396 # Loop over the parameters. 397 for param in ['S2', 'S2f', 'S2s', 'te', 'tf', 'ts', 'Rex']: 398 # The file name. 399 file_name = dir + '/' + param + '.out' 400 401 # Test if the file exists. 402 if not access(file_name, F_OK): 403 raise RelaxFileError, ('Dasha', file_name) 404 405 # Scaling. 406 if param in ['te', 'tf', 'ts']: 407 scaling = 1e-9 408 elif param == 'Rex': 409 scaling = 1.0 / (2.0 * pi * self.relax.data.frq[self.run][0]) ** 2 410 else: 411 scaling = 1.0 412 413 # Set the values. 414 self.relax.generic.value.read(self.run, param=param, scaling=scaling, file=file_name, num_col=0, name_col=None, data_col=1, error_col=2) 415 416 # Clean up of non-existant parameters (set the parameter to None!). 417 for i in xrange(len(self.relax.data.res[self.run])): 418 # Skip unselected residues. 419 if not self.relax.data.res[self.run][i].select: 420 continue 421 422 # Skip the residue (don't set the parameter to None) if the parameter exists in the model. 423 if param in self.relax.data.res[self.run][i].params: 424 continue 425 426 # Set the parameter to None. 427 setattr(self.relax.data.res[self.run][i], lower(param), None) 428 429 # Extract the chi-squared values. 430 file_name = dir + '/chi2.out' 431 432 # Test if the file exists. 433 if not access(file_name, F_OK): 434 raise RelaxFileError, ('Dasha', file_name) 435 436 # Set the values. 437 self.relax.generic.value.read(self.run, param='chi2', file=file_name, num_col=0, name_col=None, data_col=1, error_col=2)
438