generic_fns.palmer

1 ############################################################################### 2 # # 3 # Copyright (C) 2003-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 Art Palmer's Modelfree 4 program.""" 25 26 27 # Python module imports. 28 from math import pi 29 from os import F_OK, access, chdir, chmod, getcwd, listdir, popen3, remove, sep, system 30 from re import match, search 31 from string import count, find, split 32 import sys 33 34 # relax module imports. 35 from generic_fns.mol_res_spin import exists_mol_res_spin_data, spin_loop 36 from generic_fns import diffusion_tensor, pipes 37 from physical_constants import return_gyromagnetic_ratio 38 from relax_errors import RelaxDirError, RelaxFileError, RelaxNoModelError, RelaxNoPdbError, RelaxNoSequenceError 39 from relax_io import mkdir_nofail, open_write_file, test_binary 40 from specific_fns.setup import model_free_obj 41 42

43 -def __deselect_spins():

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

68 69

70 -def create(dir=None, binary=None, diff_search=None, sims=None, sim_type=None, trim=None, steps=None, heteronuc_type=None, atom1=None, atom2=None, spin_id=None, force=False, constraints=True):

71 """Create the Modelfree4 input files. 72 73 The following files are created: 74 - dir/mfin 75 - dir/mfdata 76 - dir/mfpar 77 - dir/mfmodel 78 - dir/run.sh 79 80 @keyword dir: The optional directory to place the files into. If None, then the 81 files will be placed into a directory named after the current data 82 pipe. 83 @type dir: str or None 84 @keyword binary: The name of the Modelfree4 binary file. This can include the path 85 to the binary. 86 @type binary: str 87 @keyword diff_search: The diffusion tensor search algorithm (see the Modelfree4 manual for 88 details). 89 @type diff_search: str 90 @keyword sims: The number of Monte Carlo simulations to perform. 91 @type sims: int 92 @keyword sim_type: The type of simulation to perform (see the Modelfree4 manual for 93 details). 94 @type sim_type: str 95 @keyword trim: Trimming of the Monte Carlo simulations (see the Modelfree4 manual 96 for details). 97 @type trim: int 98 @keyword steps: The grid search size (see the Modelfree4 manual for details). 99 @type steps: int 100 @keyword heteronuc_type: The Modelfree4 three letter code for the heteronucleus type, e.g. 101 '15N', '13C', etc. 102 @type heteronuc_type: str 103 @keyword atom1: The name of the heteronucleus in the PDB file. 104 @type atom1: str 105 @keyword atom2: The name of the proton in the PDB file. 106 @type atom2: str 107 @keyword spin_id: The spin identification string. 108 @type spin_id: str 109 @keyword force: A flag which if True will cause all pre-existing files to be 110 overwritten. 111 @type force: bool 112 @keyword constraints: A flag which if True will result in constrained optimisation. 113 @type constraints: bool 114 """ 115 116 # Test if the current pipe exists. 117 pipes.test() 118 119 # Test if sequence data is loaded. 120 if not exists_mol_res_spin_data(): 121 raise RelaxNoSequenceError 122 123 # Test if the PDB file is loaded (for the spheroid and ellipsoid). 124 if hasattr(cdp, 'diff_tensor') and not cdp.diff_tensor.type == 'sphere' and not hasattr(cdp, 'structure'): 125 raise RelaxNoPdbError 126 127 # Deselect certain spins. 128 __deselect_spins() 129 130 # Directory creation. 131 if dir == None: 132 dir = pipes.cdp_name() 133 mkdir_nofail(dir, verbosity=0) 134 135 # Number of field strengths and values. 136 frq = [] 137 for ri_id in cdp.ri_ids: 138 # New frequency. 139 if cdp.frq[ri_id] not in frq: 140 frq.append(cdp.frq[ri_id]) 141 142 # The 'mfin' file. 143 mfin = open_write_file('mfin', dir, force) 144 create_mfin(mfin, diff_search=diff_search, sims=sims, sim_type=sim_type, trim=trim, num_frq=len(frq), frq=frq) 145 mfin.close() 146 147 # Open the 'mfdata', 'mfmodel', and 'mfpar' files. 148 mfdata = open_write_file('mfdata', dir, force) 149 mfmodel = open_write_file('mfmodel', dir, force) 150 mfpar = open_write_file('mfpar', dir, force) 151 152 # Loop over the sequence. 153 for spin, mol_name, res_num, res_name, id in spin_loop(spin_id, full_info=True, return_id=True): 154 # Skip deselected spins. 155 if not spin.select: 156 continue 157 158 # The 'mfdata' file. 159 if not create_mfdata(mfdata, spin=spin, spin_id=id, num_frq=len(frq), frq=frq): 160 continue 161 162 # The 'mfmodel' file. 163 create_mfmodel(mfmodel, spin=spin, spin_id=id, steps=steps, constraints=constraints) 164 165 # The 'mfpar' file. 166 create_mfpar(mfpar, spin=spin, spin_id=id, res_num=res_num, atom1=atom1, atom2=atom2) 167 168 # Close the 'mfdata', 'mfmodel', and 'mfpar' files. 169 mfdata.close() 170 mfmodel.close() 171 mfpar.close() 172 173 # The 'run.sh' script. 174 run = open_write_file('run.sh', dir, force) 175 create_run(run, binary=binary, dir=dir) 176 run.close() 177 chmod(dir + sep+'run.sh', 0755)

178 179

180 -def create_mfdata(file, spin=None, spin_id=None, num_frq=None, frq=None):

181 """Create the Modelfree4 input file 'mfmodel'. 182 183 @param file: The writable file object. 184 @type file: file object 185 @param spin: The spin container. 186 @type spin: SpinContainer instance 187 @param spin_id: The spin identification string. 188 @type spin_id str 189 @keyword num_frq: The number of spectrometer frequencies relaxation data was collected at. 190 @type num_frq: int 191 @keyword frq: The spectrometer frequencies. 192 @type frq: list of float 193 @return: True if file data is written, False otherwise. 194 @rtype: bool 195 """ 196 197 # Spin title. 198 file.write("\nspin " + spin_id + "\n") 199 200 # Data written flag. 201 written = False 202 203 # Loop over the frequencies. 204 for j in xrange(num_frq): 205 # Set the data to None. 206 r1, r2, noe = None, None, None 207 208 # Loop over the relaxation data. 209 for ri_id in cdp.ri_ids: 210 # The frequency does not match. 211 if frq[j] != cdp.frq[ri_id]: 212 continue 213 214 # Find the corresponding R1. 215 if cdp.ri_type[ri_id] == 'R1': 216 r1 = spin.ri_data[ri_id] 217 r1_err = spin.ri_data_err[ri_id] 218 219 # Find the corresponding R2. 220 elif cdp.ri_type[ri_id] == 'R2': 221 r2 = spin.ri_data[ri_id] 222 r2_err = spin.ri_data_err[ri_id] 223 224 # Find the corresponding NOE. 225 elif cdp.ri_type[ri_id] == 'NOE': 226 noe = spin.ri_data[ri_id] 227 noe_err = spin.ri_data_err[ri_id] 228 229 # Test if the R1 exists for this frequency, otherwise skip the data. 230 if r1: 231 file.write('%-7s%-10.3f%20.15f%20.15f %-3i\n' % ('R1', frq[j]*1e-6, r1, r1_err, 1)) 232 else: 233 file.write('%-7s%-10.3f%20.15f%20.15f %-3i\n' % ('R1', frq[j]*1e-6, 0, 0, 0)) 234 235 # Test if the R2 exists for this frequency, otherwise skip the data. 236 if r2: 237 file.write('%-7s%-10.3f%20.15f%20.15f %-3i\n' % ('R2', frq[j]*1e-6, r2, r2_err, 1)) 238 else: 239 file.write('%-7s%-10.3f%20.15f%20.15f %-3i\n' % ('R2', frq[j]*1e-6, 0, 0, 0)) 240 241 # Test if the NOE exists for this frequency, otherwise skip the data. 242 if noe: 243 file.write('%-7s%-10.3f%20.15f%20.15f %-3i\n' % ('NOE', frq[j]*1e-6, noe, noe_err, 1)) 244 else: 245 file.write('%-7s%-10.3f%20.15f%20.15f %-3i\n' % ('NOE', frq[j]*1e-6, 0, 0, 0)) 246 247 written = True 248 249 return written

250 251

252 -def create_mfin(file, diff_search=None, sims=None, sim_type=None, trim=None, num_frq=None, frq=None):

253 """Create the Modelfree4 input file 'mfin'. 254 255 @param file: The writable file object. 256 @type file: file object 257 @keyword diff_search: The diffusion tensor search algorithm (see the Modelfree4 manual for 258 details). 259 @type diff_search: str 260 @keyword sims: The number of Monte Carlo simulations to perform. 261 @type sims: int 262 @keyword sim_type: The type of simulation to perform (see the Modelfree4 manual for 263 details). 264 @type sim_type: str 265 @keyword trim: Trimming of the Monte Carlo simulations (see the Modelfree4 manual for 266 details). 267 @type trim: int 268 @keyword num_frq: The number of spectrometer frequencies relaxation data was collected at. 269 @type num_frq: int 270 @keyword frq: The spectrometer frequencies. 271 @type frq: list of float 272 """ 273 274 # Set the diffusion tensor specific values. 275 if cdp.diff_tensor.type == 'sphere': 276 diff = 'isotropic' 277 algorithm = 'brent' 278 tm = cdp.diff_tensor.tm / 1e-9 279 dratio = 1 280 theta = 0 281 phi = 0 282 elif cdp.diff_tensor.type == 'spheroid': 283 diff = 'axial' 284 algorithm = 'powell' 285 tm = cdp.diff_tensor.tm / 1e-9 286 dratio = cdp.diff_tensor.Dratio 287 theta = cdp.diff_tensor.theta * 360.0 / (2.0 * pi) 288 phi = cdp.diff_tensor.phi * 360.0 / (2.0 * pi) 289 elif cdp.diff_tensor.type == 'ellipsoid': 290 diff = 'anisotropic' 291 algorithm = 'powell' 292 tm = cdp.diff_tensor.tm / 1e-9 293 dratio = 0 294 theta = 0 295 phi = 0 296 297 # Add the main options. 298 file.write("optimization tval\n\n") 299 file.write("seed 0\n\n") 300 file.write("search grid\n\n") 301 302 # Diffusion type. 303 if cdp.diff_tensor.fixed: 304 algorithm = 'fix' 305 306 file.write("diffusion " + diff + " " + diff_search + "\n\n") 307 file.write("algorithm " + algorithm + "\n\n") 308 309 # Monte Carlo simulations. 310 if sims: 311 file.write("simulations " + sim_type + " " + repr(sims) + " " + repr(trim) + "\n\n") 312 else: 313 file.write("simulations none\n\n") 314 315 selection = 'none' # To be changed. 316 file.write("selection " + selection + "\n\n") 317 file.write("sim_algorithm " + algorithm + "\n\n") 318 319 file.write("fields " + repr(num_frq)) 320 for val in frq: 321 file.write(" " + repr(val*1e-6)) 322 file.write("\n") 323 324 # tm. 325 file.write('%-7s' % 'tm') 326 file.write('%14.3f' % tm) 327 file.write('%2i' % 1) 328 file.write('%3i' % 0) 329 file.write('%5i' % 5) 330 file.write('%6i' % 15) 331 file.write('%4i\n' % 20) 332 333 # dratio. 334 file.write('%-7s' % 'Dratio') 335 file.write('%14s' % dratio) 336 file.write('%2i' % 1) 337 file.write('%3i' % 0) 338 file.write('%5i' % 0) 339 file.write('%6i' % 2) 340 file.write('%4i\n' % 5) 341 342 # theta. 343 file.write('%-7s' % 'Theta') 344 file.write('%14s' % theta) 345 file.write('%2i' % 1) 346 file.write('%3i' % 0) 347 file.write('%5i' % 0) 348 file.write('%6i' % 180) 349 file.write('%4i\n' % 10) 350 351 # phi. 352 file.write('%-7s' % 'Phi') 353 file.write('%14s' % phi) 354 file.write('%2i' % 1) 355 file.write('%3i' % 0) 356 file.write('%5i' % 0) 357 file.write('%6i' % 360) 358 file.write('%4i\n' % 10)

359 360

361 -def create_mfmodel(file, spin=None, spin_id=None, steps=None, constraints=None):

362 """Create the Modelfree4 input file 'mfmodel'. 363 364 @param file: The writable file object. 365 @type file: file object 366 @keyword spin: The spin container. 367 @type spin: SpinContainer instance 368 @keyword spin_id: The spin identification string. 369 @type spin_id str 370 @keyword steps: The grid search size (see the Modelfree4 manual for details). 371 @type steps: int 372 @keyword constraints: A flag which if True will result in constrained optimisation. 373 @type constraints: bool 374 """ 375 376 # Spin title. 377 file.write("\nspin " + spin_id + "\n") 378 379 # tloc. 380 file.write('%-3s%-6s%-6.1f' % ('M1', 'tloc', 0)) 381 if 'tm' in spin.params: 382 file.write('%-4i' % 1) 383 else: 384 file.write('%-4i' % 0) 385 386 if constraints: 387 file.write('%-2i' % 2) 388 else: 389 file.write('%-2i' % 0) 390 391 file.write('%11.3f%12.3f %-4s\n' % (0, 20, steps)) 392 393 # Theta. 394 file.write('%-3s%-6s%-6.1f' % ('M1', 'Theta', 0)) 395 file.write('%-4i' % 0) 396 397 if constraints: 398 file.write('%-2i' % 2) 399 else: 400 file.write('%-2i' % 0) 401 402 file.write('%11.3f%12.3f %-4s\n' % (0, 90, steps)) 403 404 # S2f. 405 file.write('%-3s%-6s%-6.1f' % ('M1', 'Sf2', 1)) 406 if 's2f' in spin.params: 407 file.write('%-4i' % 1) 408 else: 409 file.write('%-4i' % 0) 410 411 if constraints: 412 file.write('%-2i' % 2) 413 else: 414 file.write('%-2i' % 0) 415 416 file.write('%11.3f%12.3f %-4s\n' % (0, 1, steps)) 417 418 # S2s. 419 file.write('%-3s%-6s%-6.1f' % ('M1', 'Ss2', 1)) 420 if 's2s' in spin.params or 's2' in spin.params: 421 file.write('%-4i' % 1) 422 else: 423 file.write('%-4i' % 0) 424 425 if constraints: 426 file.write('%-2i' % 2) 427 else: 428 file.write('%-2i' % 0) 429 430 file.write('%11.3f%12.3f %-4s\n' % (0, 1, steps)) 431 432 # te. 433 file.write('%-3s%-6s%-6.1f' % ('M1', 'te', 0)) 434 if 'te' in spin.params or 'ts' in spin.params: 435 file.write('%-4i' % 1) 436 else: 437 file.write('%-4i' % 0) 438 439 if constraints: 440 file.write('%-2i' % 2) 441 else: 442 file.write('%-2i' % 0) 443 444 file.write('%11.3f%12.3f %-4s\n' % (0, 10000, steps)) 445 446 # Rex. 447 file.write('%-3s%-6s%-6.1f' % ('M1', 'rex', 0)) 448 if 'rex' in spin.params: 449 file.write('%-4i' % 1) 450 else: 451 file.write('%-4i' % 0) 452 453 if constraints: 454 file.write('%-2i' % -1) 455 else: 456 file.write('%-2i' % 0) 457 458 file.write('%11.3f%12.3f %-4s\n' % (0, 20, steps))

459 460

461 -def create_mfpar(file, spin=None, spin_id=None, res_num=None, atom1=None, atom2=None):

462 """Create the Modelfree4 input file 'mfpar'. 463 464 @param file: The writable file object. 465 @type file: file object 466 @keyword spin: The spin container. 467 @type spin: SpinContainer instance 468 @keyword spin_id: The spin identification string. 469 @type spin_id str 470 @keyword res_num: The residue number from the PDB file corresponding to the spin. 471 @type res_num: int 472 @keyword atom1: The name of the heteronucleus in the PDB file. 473 @type atom1: str 474 @keyword atom2: The name of the proton in the PDB file. 475 @type atom2: str 476 """ 477 478 # Spin title. 479 file.write("\nspin " + spin_id + "\n") 480 481 file.write('%-14s' % "constants") 482 file.write('%-6i' % res_num) 483 file.write('%-7s' % spin.heteronuc_type) 484 file.write('%-8.4f' % (return_gyromagnetic_ratio(spin.heteronuc_type) / 1e7)) 485 file.write('%-8.3f' % (spin.r * 1e10)) 486 file.write('%-8.3f\n' % (spin.csa * 1e6)) 487 488 file.write('%-10s' % "vector") 489 file.write('%-4s' % atom1) 490 file.write('%-4s\n' % atom2)

491 492

493 -def create_run(file, binary=None, dir=None):

494 """Create the script 'run.sh' for the execution of Modelfree4. 495 496 @param file: The writable file object. 497 @type file: file object 498 @keyword binary: The name of the Modelfree4 binary file. This can include the path to the 499 binary. 500 @type binary: str 501 @keyword dir: The directory to copy the PDB file to. 502 @type dir: str 503 """ 504 505 file.write("#! /bin/sh\n") 506 file.write(binary + " -i mfin -d mfdata -p mfpar -m mfmodel -o mfout -e out") 507 if cdp.diff_tensor.type != 'sphere': 508 # Copy the first pdb file to the model directory so there are no problems with the existance of *.rotate files. 509 mol = cdp.structure.structural_data[0].mol[0] 510 system('cp ' + mol.file_path + sep + mol.file_name + ' ' + dir) 511 file.write(" -s " + mol.file_name) 512 file.write("\n")

513 514

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

516 """Execute Modelfree4. 517 518 BUG: Control-C during execution causes the cwd to stay as dir. 519 520 521 @param dir: The optional directory where the script is located. 522 @type dir: str or None 523 @param force: A flag which if True will cause any pre-existing files to be overwritten by 524 Modelfree4. 525 @type force: bool 526 @param binary: The name of the Modelfree4 binary file. This can include the path to the 527 binary. 528 @type binary: str 529 """ 530 531 # The current directory. 532 orig_dir = getcwd() 533 534 # The directory. 535 if dir == None: 536 dir = pipes.cdp_name() 537 if not access(dir, F_OK): 538 raise RelaxDirError('Modelfree4', dir) 539 540 # Change to this directory. 541 chdir(dir) 542 543 # Catch failures and return to the correct directory. 544 try: 545 # Test if the 'mfin' input file exists. 546 if not access('mfin', F_OK): 547 raise RelaxFileError('mfin input', 'mfin') 548 549 # Test if the 'mfdata' input file exists. 550 if not access('mfdata', F_OK): 551 raise RelaxFileError('mfdata input', 'mfdata') 552 553 # Test if the 'mfmodel' input file exists. 554 if not access('mfmodel', F_OK): 555 raise RelaxFileError('mfmodel input', 'mfmodel') 556 557 # Test if the 'mfpar' input file exists. 558 if not access('mfpar', F_OK): 559 raise RelaxFileError('mfpar input', 'mfpar') 560 561 # Test if the 'PDB' input file exists. 562 if cdp.diff_tensor.type != 'sphere': 563 pdb = cdp.structure.structural_data[0].mol[0].file_name 564 if not access(pdb, F_OK): 565 raise RelaxFileError('PDB', pdb) 566 else: 567 pdb = None 568 569 # Remove the file 'mfout' and '*.out' if the force flag is set. 570 if force: 571 for file in listdir(getcwd()): 572 if search('out$', file) or search('rotate$', file): 573 remove(file) 574 575 # Test the binary file string corresponds to a valid executable. 576 test_binary(binary) 577 578 # Execute Modelfree4. 579 if pdb: 580 cmd = binary + ' -i mfin -d mfdata -p mfpar -m mfmodel -o mfout -e out -s ' + pdb 581 else: 582 cmd = binary + ' -i mfin -d mfdata -p mfpar -m mfmodel -o mfout -e out' 583 stdin, stdout, stderr = popen3(cmd) 584 585 # Close the pipe. 586 stdin.close() 587 588 # Write to stdout and stderr. 589 for line in stdout.readlines(): 590 sys.stdout.write(line) 591 for line in stderr.readlines(): 592 sys.stderr.write(line) 593 594 # Failure. 595 except: 596 # Change back to the original directory. 597 chdir(orig_dir) 598 599 # Reraise the error. 600 raise 601 602 # Change back to the original directory. 603 chdir(orig_dir)

604 605

606 -def extract(dir, spin_id=None):

607 """Extract the Modelfree4 results out of the 'mfout' file. 608 609 @param dir: The directory containing the 'mfout' file. 610 @type dir: str or None 611 @keyword spin_id: The spin identification string. 612 @type spin_id: str or None 613 """ 614 615 # Test if sequence data is loaded. 616 if not exists_mol_res_spin_data(): 617 raise RelaxNoSequenceError 618 619 # The directory. 620 if dir == None: 621 dir = pipes.cdp_name() 622 if not access(dir, F_OK): 623 raise RelaxDirError('Modelfree4', dir) 624 625 # Test if the file exists. 626 if not access(dir + sep+'mfout', F_OK): 627 raise RelaxFileError('Modelfree4', dir + sep+'mfout') 628 629 # Determine the parameter set. 630 model_type = model_free_obj._determine_model_type() 631 632 # Open the file. 633 mfout_file = open(dir + sep+'mfout', 'r') 634 mfout_lines = mfout_file.readlines() 635 mfout_file.close() 636 637 # Get the section line positions of the mfout file. 638 global_chi2_pos, diff_pos, s2_pos, s2f_pos, s2s_pos, te_pos, rex_pos, chi2_pos = line_positions(mfout_lines) 639 640 # Find out if simulations were carried out. 641 sims = 0 642 for i in xrange(len(mfout_lines)): 643 if search('_iterations', mfout_lines[i]): 644 row = split(mfout_lines[i]) 645 sims = int(row[1]) 646 647 # Global data. 648 if model_type in ['all', 'diff']: 649 # Global chi-squared. 650 row = split(mfout_lines[global_chi2_pos]) 651 cdp.chi2 = float(row[1]) 652 653 # Spherical diffusion tensor. 654 if cdp.diff_tensor.type == 'sphere': 655 # Split the lines. 656 tm_row = split(mfout_lines[diff_pos]) 657 658 # Set the params. 659 cdp.diff_tensor.tm = float(tm_row[2]) 660 661 # Spheroid diffusion tensor. 662 else: 663 # Split the lines. 664 tm_row = split(mfout_lines[diff_pos]) 665 dratio_row = split(mfout_lines[diff_pos+1]) 666 theta_row = split(mfout_lines[diff_pos+2]) 667 phi_row = split(mfout_lines[diff_pos+3]) 668 669 # Set the params. 670 diffusion_tensor.set([float(tm_row[2]), float(dratio_row[2]), float(theta_row[2])*2.0*pi/360.0, float(phi_row[2])*2.0*pi/360.0], ['tm', 'Dratio', 'theta', 'phi']) 671 672 # Loop over the sequence. 673 pos = 0 674 for spin, mol_name, res_num, res_name in spin_loop(spin_id, full_info=True): 675 # Skip deselected residues. 676 if not spin.select: 677 continue 678 679 # Get the residue number from the mfout file. 680 mfout_res_num = int(split(mfout_lines[s2_pos + pos])[0]) 681 682 # Skip the spin if the residue doesn't match. 683 if mfout_res_num != res_num: 684 continue 685 686 # Test that the model has been set (needed to differentiate between te and ts). 687 if not hasattr(spin, 'model'): 688 raise RelaxNoModelError 689 690 # Get the S2 data. 691 if 's2' in spin.params: 692 spin.s2, spin.s2_err = get_mf_data(mfout_lines, s2_pos + pos) 693 694 # Get the S2f data. 695 if 's2f' in spin.params or 's2s' in spin.params: 696 spin.s2f, spin.s2f_err = get_mf_data(mfout_lines, s2f_pos + pos) 697 698 # Get the S2s data. 699 if 's2f' in spin.params or 's2s' in spin.params: 700 spin.s2s, spin.s2s_err = get_mf_data(mfout_lines, s2s_pos + pos) 701 702 # Get the te data. 703 if 'te' in spin.params: 704 spin.te, spin.te_err = get_mf_data(mfout_lines, te_pos + pos) 705 spin.te = spin.te / 1e12 706 spin.te_err = spin.te_err / 1e12 707 708 # Get the ts data. 709 if 'ts' in spin.params: 710 spin.ts, spin.ts_err = get_mf_data(mfout_lines, te_pos + pos) 711 spin.ts = spin.ts / 1e12 712 spin.ts_err = spin.ts_err / 1e12 713 714 # Get the Rex data. 715 if 'rex' in spin.params: 716 spin.rex, spin.rex_err = get_mf_data(mfout_lines, rex_pos + pos) 717 spin.rex = spin.rex / (2.0 * pi * cdp.frq[cdp.ri_ids[0]])**2 718 spin.rex_err = spin.rex_err / (2.0 * pi * cdp.frq[cdp.ri_ids[0]])**2 719 720 # Get the chi-squared data. 721 if not sims: 722 row = split(mfout_lines[chi2_pos + pos]) 723 spin.chi2 = float(row[1]) 724 else: 725 # The mfout chi2 position (with no sims) plus 2 (for the extra XML) plus the residue position times 22 (because of the simulated SSE rows). 726 row = split(mfout_lines[chi2_pos + 2 + 22*pos]) 727 spin.chi2 = float(row[1]) 728 729 # Increment the residue position. 730 pos = pos + 1

731 732

733 -def get_mf_data(mfout_lines, pos):

734 """Extract the model-free data from the given position of the mfout file. 735 736 This method is designed to catch a number of bugs in Modelfree4's mfout file. 737 738 The first bug is the presence of a series of '*' characters causing a fusion of two columns. 739 This is handled by splitting by the '*' char and then returning the first element. 740 741 The second bug is when the floating point number is too big to fit into Modelfree4's string 742 format limit of 15.3f. This results in a results line such as: 743 744 246 10000.00019682363392.000 1 0.000 0.000 0.000 0.000 745 746 This is caught by scanning for two '.' characters in the column, and handled by assuming 747 that every floating point number will have three decimal characters. 748 749 @param mfout_lines: A list of all the lines of the mfout file. 750 @type mfout_lines: list of str 751 @param pos: The mfout line position. 752 @type pos: int 753 @return: The value and error. 754 @rtype: tuple of 2 floats 755 """ 756 757 # Split the line up. 758 row = split(mfout_lines[pos]) 759 760 # The value and error, assuming a bug free mfout file. 761 val = row[1] 762 err = row[4] 763 764 # The Modelfree4 '*' column fusion bug. 765 if search('\*', val) or search('\*', err): 766 # Split by the '*' character. 767 val_row = split(val, '*') 768 err_row = split(err, '*') 769 770 # The value and error (the first elements). 771 val = val_row[0] 772 err = err_row[0] 773 774 # The Modelfree4 large float column fusion bug. 775 new_row = [] 776 fused = False 777 for element in row: 778 # Count the number of '.' characters. 779 num = count(element, '.') 780 781 # Catch two or more '.' characters. 782 if num > 1: 783 # Set the fused flag. 784 fused = True 785 786 # Loop over each fused number. 787 for i in xrange(num): 788 # Find the index of the first '.'. 789 index = find(element, '.') 790 791 # The first number (index + decimal point + 3 decimal chars). 792 new_row.append(element[0:index+4]) 793 794 # Strip the first number from the element for the next loop iteration. 795 element = element[index+4:] 796 797 # Otherwise the column element is fine. 798 else: 799 new_row.append(element) 800 801 # Bug has been caught. 802 if fused: 803 val = new_row[1] 804 err = new_row[4] 805 806 # Return the value and error, as floats. 807 return float(val), float(err)

808 809

810 -def line_positions(mfout_lines):

811 """Function for getting the section positions (line number) of the mfout file. 812 813 @param mfout_lines: A list of all the lines of the mfout file. 814 @type mfout_lines: list of str 815 @return: The line indices where the s2, s2f, s2s, te, rex, and chi2 sections 816 start in the mfout file. 817 @rtype: tuple of int 818 """ 819 820 # Loop over the file. 821 i = 0 822 while i < len(mfout_lines): 823 # Global chi2. 824 if match('data_chi_square', mfout_lines[i]): 825 global_chi2_pos = i + 1 826 827 # Diffusion tensor. 828 if match('data_diffusion_tensor', mfout_lines[i]): 829 diff_pos = i + 3 830 831 # Model-free data. 832 if match('data_model_1', mfout_lines[i]): 833 # Shift down two lines (to avoid the lines not starting with a space). 834 i = i + 2 835 836 # Walk through all the data. 837 while True: 838 # Break once the end of the data section is reached. 839 if not mfout_lines[i] == '\n' and not search('^ ', mfout_lines[i]): 840 break 841 842 # Split the line up. 843 row = split(mfout_lines[i]) 844 845 # S2 position (skip the heading and move to the first residue). 846 if len(row) == 2 and row[0] == 'S2': 847 s2_pos = i + 1 848 849 # S2f position (skip the heading and move to the first residue). 850 if len(row) == 2 and row[0] == 'S2f': 851 s2f_pos = i + 1 852 853 # S2s position (skip the heading and move to the first residue). 854 if len(row) == 2 and row[0] == 'S2s': 855 s2s_pos = i + 1 856 857 # te position (skip the heading and move to the first residue). 858 if len(row) == 2 and row[0] == 'te': 859 te_pos = i + 1 860 861 # Rex position (skip the heading and move to the first residue). 862 if len(row) == 2 and row[0] == 'Rex': 863 rex_pos = i + 1 864 865 # Move to the next line number. 866 i = i + 1 867 868 # Chi-squared values. 869 if match('data_sse', mfout_lines[i]): 870 chi2_pos = i + 3 871 872 # Move to the next line number. 873 i = i + 1 874 875 # Return the positions. 876 return global_chi2_pos, diff_pos, s2_pos, s2f_pos, s2s_pos, te_pos, rex_pos, chi2_pos

877

Source Code for Module generic_fns.palmer