generic_fns.dasha

1 ############################################################################### 2 # # 3 # Copyright (C) 2005-2012 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 # Module docstring. 24 """Module for interfacing with Dasha.""" 25 26 # Python module imports. 27 from math import pi 28 from os import F_OK, access, chdir, getcwd, popen3, sep 29 from string import lower 30 import sys 31 32 # relax module imports. 33 from generic_fns import angles, diffusion_tensor, pipes, relax_data, value 34 from generic_fns.mol_res_spin import exists_mol_res_spin_data, first_residue_num, last_residue_num, residue_loop, spin_loop 35 from relax_errors import RelaxDirError, RelaxError, RelaxFileError, RelaxNoPdbError, RelaxNoSequenceError, RelaxNoTensorError 36 from relax_io import mkdir_nofail, open_write_file, test_binary 37 from specific_fns.setup import model_free_obj 38 39

40 -def __deselect_spins():

41 """Deselect spins with no or too little data, that are overfitting, etc.""" 42 43 # Test if sequence data exists. 44 if not exists_mol_res_spin_data(): 45 raise RelaxNoSequenceError 46 47 # Is structural data required? 48 need_vect = False 49 if hasattr(cdp, 'diff_tensor') and (cdp.diff_tensor.type == 'spheroid' or cdp.diff_tensor.type == 'ellipsoid'): 50 need_vect = True 51 52 # Loop over the sequence. 53 for spin in spin_loop(): 54 # Relaxation data must exist! 55 if not hasattr(spin, 'ri_data'): 56 spin.select = False 57 58 # Require 3 or more relaxation data points. 59 elif len(spin.ri_data) < 3: 60 spin.select = False 61 62 # Require at least as many data points as params to prevent over-fitting. 63 elif hasattr(spin, 'params') and spin.params and len(spin.params) > len(spin.ri_data): 64 spin.select = False

65 66

67 -def create(algor='LM', dir=None, force=False):

68 """Create the Dasha script file 'dasha_script' for controlling the program. 69 70 @keyword algor: The optimisation algorithm to use. This can be the Levenberg-Marquardt algorithm 'LM' or the Newton-Raphson algorithm 'NR'. 71 @type algor: str 72 @keyword dir: The optional directory to place the script into. 73 @type dir: str or None 74 @keyword force: A flag which if True will cause any pre-existing file to be overwritten. 75 @type force: bool 76 """ 77 78 # Test if the current pipe exists. 79 pipes.test() 80 81 # Test if sequence data is loaded. 82 if not exists_mol_res_spin_data(): 83 raise RelaxNoSequenceError 84 85 # Determine the parameter set. 86 model_type = model_free_obj._determine_model_type() 87 88 # Test if diffusion tensor data for the data_pipe exists. 89 if model_type != 'local_tm' and not hasattr(cdp, 'diff_tensor'): 90 raise RelaxNoTensorError('diffusion') 91 92 # Test if the PDB file has been loaded (for the spheroid and ellipsoid). 93 if model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere' and not hasattr(cdp, 'structure'): 94 raise RelaxNoPdbError 95 96 # Test the optimisation algorithm. 97 if algor not in ['LM', 'NR']: 98 raise RelaxError("The Dasha optimisation algorithm '%s' is unknown, it should either be 'LM' or 'NR'." % algor) 99 100 # Multiple spins per residue not allowed. 101 for residue in residue_loop(): 102 # Test the number of spins. 103 if len(residue.spin) > 1: 104 raise RelaxError("More than one spin per residue is not supported.") 105 106 # Deselect certain spins. 107 __deselect_spins() 108 109 # Directory creation. 110 if dir == None: 111 dir = pipes.cdp_name() 112 mkdir_nofail(dir, verbosity=0) 113 114 # Calculate the angle alpha of the XH vector in the spheroid diffusion frame. 115 if cdp.diff_tensor.type == 'spheroid': 116 angles.spheroid_frame() 117 118 # Calculate the angles theta and phi of the XH vector in the ellipsoid diffusion frame. 119 elif cdp.diff_tensor.type == 'ellipsoid': 120 angles.ellipsoid_frame() 121 122 # The 'dasha_script' file. 123 script = open_write_file(file_name='dasha_script', dir=dir, force=force) 124 create_script(script, model_type, algor) 125 script.close()

126 127

128 -def create_script(file, model_type, algor):

129 """Create the Dasha script file. 130 131 @param file: The opened file descriptor. 132 @type file: file object 133 @param model_type: The model-free model type. 134 @type model_type: str 135 @param algor: The optimisation algorithm to use. This can be the Levenberg-Marquardt algorithm 'LM' or the Newton-Raphson algorithm 'NR'. 136 @type algor: str 137 """ 138 139 # Delete all data. 140 file.write("# Delete all data.\n") 141 file.write("del 1 10000\n") 142 143 # Nucleus type. 144 file.write("\n# Nucleus type.\n") 145 nucleus = None 146 for spin in spin_loop(): 147 # Can only handle one spin type. 148 if nucleus and spin.heteronuc_type != nucleus: 149 raise RelaxError("The nuclei '%s' and '%s' do not match, relax can only handle one nucleus type in Dasha." % (nucleus, spin.heteronuc_type)) 150 151 # Set the nucleus. 152 if not nucleus: 153 nucleus = spin.heteronuc_type 154 155 # Convert the name and write it. 156 if nucleus == '15N': 157 nucleus = 'N15' 158 elif nucleus == '13C': 159 nucleus = 'C13' 160 else: 161 raise RelaxError("Cannot handle the nucleus type '%s' within Dasha." % nucleus) 162 file.write("set nucl %s\n" % nucleus) 163 164 # Number of frequencies. 165 file.write("\n# Number of frequencies.\n") 166 file.write("set n_freq %s\n" % relax_data.num_frq()) 167 168 # Frequency values. 169 file.write("\n# Frequency values.\n") 170 count = 1 171 for frq in relax_data.frq_loop(): 172 file.write("set H1_freq %s %s\n" % (frq / 1e6, count)) 173 count += 1 174 175 # Set the diffusion tensor. 176 file.write("\n# Set the diffusion tensor.\n") 177 if model_type != 'local_tm': 178 # Sphere. 179 if cdp.diff_tensor.type == 'sphere': 180 file.write("set tr %s\n" % (cdp.diff_tensor.tm / 1e-9)) 181 182 # Spheroid. 183 elif cdp.diff_tensor.type == 'spheroid': 184 file.write('set tr %s\n' % (cdp.diff_tensor.tm / 1e-9)) 185 186 # Ellipsoid. 187 elif cdp.diff_tensor.type == 'ellipsoid': 188 # Get the eigenvales. 189 Dx, Dy, Dz = diffusion_tensor.return_eigenvalues() 190 191 # Geometric parameters. 192 file.write("set tr %s\n" % (cdp.diff_tensor.tm / 1e-9)) 193 file.write("set D1/D3 %s\n" % (Dx / Dz)) 194 file.write("set D2/D3 %s\n" % (Dy / Dz)) 195 196 # Orientational parameters. 197 file.write("set alfa %s\n" % (cdp.diff_tensor.alpha / (2.0 * pi) * 360.0)) 198 file.write("set betta %s\n" % (cdp.diff_tensor.beta / (2.0 * pi) * 360.0)) 199 file.write("set gamma %s\n" % (cdp.diff_tensor.gamma / (2.0 * pi) * 360.0)) 200 201 # Reading the relaxation data. 202 file.write("\n# Reading the relaxation data.\n") 203 file.write("echo Reading the relaxation data.\n") 204 noe_index = 1 205 r1_index = 1 206 r2_index = 1 207 for ri_id in cdp.ri_ids: 208 # NOE. 209 if cdp.ri_type[ri_id] == 'NOE': 210 # Data set number. 211 number = noe_index 212 213 # Data type. 214 data_type = 'noe' 215 216 # Increment the data set index. 217 noe_index = noe_index + 1 218 219 # R1. 220 elif cdp.ri_type[ri_id] == 'R1': 221 # Data set number. 222 number = r1_index 223 224 # Data type. 225 data_type = '1/T1' 226 227 # Increment the data set index. 228 r1_index = r1_index + 1 229 230 # R2. 231 elif cdp.ri_type[ri_id] == 'R2': 232 # Data set number. 233 number = r2_index 234 235 # Data type. 236 data_type = '1/T2' 237 238 # Increment the data set index. 239 r2_index = r2_index + 1 240 241 # Set the data type. 242 if number == 1: 243 file.write("\nread < %s\n" % data_type) 244 else: 245 file.write("\nread < %s %s\n" % (data_type, number)) 246 247 # The relaxation data. 248 for residue in residue_loop(): 249 # Alias the spin. 250 spin = residue.spin[0] 251 252 # Skip deselected spins. 253 if not spin.select: 254 continue 255 256 # Skip and deselect spins for which relaxation data is missing. 257 if len(spin.ri_data) != len(cdp.ri_ids) or spin.ri_data[ri_id] == None: 258 spin.select = False 259 continue 260 261 # Data and errors. 262 file.write("%s %s %s\n" % (residue.num, spin.ri_data[ri_id], spin.ri_data_err[ri_id])) 263 264 # Terminate the reading. 265 file.write("exit\n") 266 267 # Individual residue optimisation. 268 if model_type == 'mf': 269 # Loop over the residues. 270 for residue in residue_loop(): 271 # Alias the spin. 272 spin = residue.spin[0] 273 274 # Skip deselected spins. 275 if not spin.select: 276 continue 277 278 # Comment. 279 file.write("\n\n\n# Residue %s\n\n" % residue.num) 280 281 # Echo. 282 file.write("echo Optimisation of residue %s\n" % residue.num) 283 284 # Select the spin. 285 file.write("\n# Select the residue.\n") 286 file.write("set cres %s\n" % residue.num) 287 288 # The angle alpha of the XH vector in the spheroid diffusion frame. 289 if cdp.diff_tensor.type == 'spheroid': 290 file.write("set teta %s\n" % spin.alpha) 291 292 # The angles theta and phi of the XH vector in the ellipsoid diffusion frame. 293 elif cdp.diff_tensor.type == 'ellipsoid': 294 file.write("\n# Setting the spherical angles of the XH vector in the ellipsoid diffusion frame.\n") 295 file.write("set teta %s\n" % spin.theta) 296 file.write("set fi %s\n" % spin.phi) 297 298 # The 'jmode'. 299 if 'ts' in spin.params: 300 jmode = 3 301 elif 'te' in spin.params: 302 jmode = 2 303 elif 's2' in spin.params: 304 jmode = 1 305 306 # Chemical exchange. 307 if 'rex' in spin.params: 308 exch = True 309 else: 310 exch = False 311 312 # Anisotropic diffusion. 313 if cdp.diff_tensor.type == 'sphere': 314 anis = False 315 else: 316 anis = True 317 318 # Axial symmetry. 319 if cdp.diff_tensor.type == 'spheroid': 320 sym = True 321 else: 322 sym = False 323 324 # Set the jmode. 325 file.write("\n# Set the jmode.\n") 326 file.write("set def jmode %s" % jmode) 327 if exch: 328 file.write(" exch") 329 if anis: 330 file.write(" anis") 331 if sym: 332 file.write(" sym") 333 file.write("\n") 334 335 # Parameter default values. 336 file.write("\n# Parameter default values.\n") 337 file.write("reset jmode %s\n" % residue.num) 338 339 # Bond length. 340 file.write("\n# Bond length.\n") 341 file.write("set r_hx %s\n" % (spin.r / 1e-10)) 342 343 # CSA value. 344 file.write("\n# CSA value.\n") 345 file.write("set csa %s\n" % (spin.csa / 1e-6)) 346 347 # Fix the tf parameter if it isn't in the model. 348 if not 'tf' in spin.params and jmode == 3: 349 file.write("\n# Fix the tf parameter.\n") 350 file.write("fix tf 0\n") 351 352 # Optimisation of all residues. 353 file.write("\n\n\n# Optimisation of all residues.\n") 354 if algor == 'LM': 355 file.write("lmin %s %s" % (first_residue_num(), last_residue_num())) 356 elif algor == 'NR': 357 file.write("min %s %s" % (first_residue_num(), last_residue_num())) 358 359 # Show the results. 360 file.write("\n# Show the results.\n") 361 file.write("echo\n") 362 file.write("show all\n") 363 364 # Write the results. 365 file.write("\n# Write the results.\n") 366 file.write("write s2.out S\n") 367 file.write("write s2f.out Sf\n") 368 file.write("write s2s.out Ss\n") 369 file.write("write te.out te\n") 370 file.write("write tf.out tf\n") 371 file.write("write ts.out ts\n") 372 file.write("write rex.out rex\n") 373 file.write("write chi2.out F\n") 374 375 else: 376 raise RelaxError("Optimisation of the parameter set '%s' currently not supported." % model_type)

377 378

379 -def execute(dir, force, binary):

380 """Execute Dasha. 381 382 @param dir: The optional directory where the script is located. 383 @type dir: str or None 384 @param force: A flag which if True will cause any pre-existing files to be overwritten by Dasha. 385 @type force: bool 386 @param binary: The name of the Dasha binary file. This can include the path to the binary. 387 @type binary: str 388 """ 389 390 # Test the binary file string corresponds to a valid executable. 391 test_binary(binary) 392 393 # The current directory. 394 orig_dir = getcwd() 395 396 # The directory. 397 if dir == None: 398 dir = pipes.cdp_name() 399 if not access(dir, F_OK): 400 raise RelaxDirError('Dasha', dir) 401 402 # Change to this directory. 403 chdir(dir) 404 405 # Catch failures and return to the correct directory. 406 try: 407 # Test if the 'dasha_script' script file exists. 408 if not access('dasha_script', F_OK): 409 raise RelaxFileError('dasha script', 'dasha_script') 410 411 # Execute Dasha. 412 stdin, stdout, stderr = popen3(binary) 413 414 # Get the contents of the script and pump it into Dasha. 415 script = open('dasha_script') 416 lines = script.readlines() 417 script.close() 418 for line in lines: 419 stdin.write(line) 420 421 # Close the pipe. 422 stdin.close() 423 424 # Write to stdout and stderr. 425 for line in stdout.readlines(): 426 sys.stdout.write(line) 427 for line in stderr.readlines(): 428 sys.stderr.write(line) 429 430 # Failure. 431 except: 432 # Change back to the original directory. 433 chdir(orig_dir) 434 435 # Reraise the error. 436 raise 437 438 # Change back to the original directory. 439 chdir(orig_dir) 440 441 # Print some blank lines (aesthetics) 442 sys.stdout.write("\n\n")

443 444

445 -def extract(dir):

446 """Extract the data from the Dasha results files. 447 448 @param dir: The optional directory where the results file is located. 449 @type dir: str or None 450 """ 451 452 # Test if sequence data is loaded. 453 if not exists_mol_res_spin_data(): 454 raise RelaxNoSequenceError 455 456 # The directory. 457 if dir == None: 458 dir = pipes.cdp_name() 459 if not access(dir, F_OK): 460 raise RelaxDirError('Dasha', dir) 461 462 # Loop over the parameters. 463 for param in ['s2', 's2f', 's2s', 'te', 'tf', 'ts', 'rex']: 464 # The file name. 465 file_name = dir + sep + param + '.out' 466 467 # Test if the file exists. 468 if not access(file_name, F_OK): 469 raise RelaxFileError('Dasha', file_name) 470 471 # Scaling. 472 if param in ['te', 'tf', 'ts']: 473 scaling = 1e-9 474 elif param == 'rex': 475 scaling = 1.0 / (2.0 * pi * cdp.frq[cdp.ri_ids[0]]) ** 2 476 else: 477 scaling = 1.0 478 479 # Set the values. 480 value.read(param=param, scaling=scaling, file=file_name, res_num_col=1, res_name_col=None, data_col=2, error_col=3) 481 482 # Clean up of non-existent parameters (set the parameter to None!). 483 for spin in spin_loop(): 484 # Skip the spin (don't set the parameter to None) if the parameter exists in the model. 485 if param in spin.params: 486 continue 487 488 # Set the parameter to None. 489 setattr(spin, lower(param), None) 490 491 # Extract the chi-squared values. 492 file_name = dir + sep+'chi2.out' 493 494 # Test if the file exists. 495 if not access(file_name, F_OK): 496 raise RelaxFileError('Dasha', file_name) 497 498 # Set the values. 499 value.read(param='chi2', file=file_name, res_num_col=1, res_name_col=None, data_col=2, error_col=3)

500

Source Code for Module generic_fns.dasha