pipe_control.palmer

1 ############################################################################### 2 # # 3 # Copyright (C) 2003-2014 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 # Module docstring. 23 """Module for interfacing with Art Palmer's Modelfree 4 program.""" 24 25 26 # Dependencies. 27 import dep_check 28 29 # Python module imports. 30 from math import pi 31 from os import F_OK, access, chdir, chmod, getcwd, listdir, remove, sep, system 32 from re import match, search 33 from stat import S_IRWXU, S_IRGRP, S_IROTH 34 PIPE, Popen = None, None 35 if dep_check.subprocess_module: 36 from subprocess import PIPE, Popen 37 import sys 38 39 # relax module imports. 40 from lib.errors import RelaxError, RelaxDirError, RelaxFileError, RelaxNoInteratomError, RelaxNoModelError, RelaxNoPdbError, RelaxNoSequenceError, RelaxNoTensorError 41 from lib.io import mkdir_nofail, open_write_file, test_binary 42 from lib.periodic_table import periodic_table 43 from pipe_control import diffusion_tensor, pipes 44 from pipe_control.interatomic import return_interatom_list 45 from pipe_control.mol_res_spin import exists_mol_res_spin_data, spin_loop 46 from pipe_control.pipes import check_pipe 47 from specific_analyses.model_free.model import determine_model_type 48 49

50 -def __deselect_spins():

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

75 76

77 -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):

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

177 178

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

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

249 250

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

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

362 363

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

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

462 463

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

465 """Create the Modelfree4 input file 'mfpar'. 466 467 @param file: The writable file object. 468 @type file: file object 469 @keyword spin: The spin container. 470 @type spin: SpinContainer instance 471 @keyword spin_id: The spin identification string. 472 @type spin_id str 473 @keyword res_num: The residue number from the PDB file corresponding to the spin. 474 @type res_num: int 475 @keyword atom1: The name of the heteronucleus in the PDB file. 476 @type atom1: str 477 @keyword atom2: The name of the proton in the PDB file. 478 @type atom2: str 479 """ 480 481 # Get the interatomic data containers. 482 interatoms = return_interatom_list(spin_id) 483 if len(interatoms) == 0: 484 raise RelaxNoInteratomError 485 elif len(interatoms) > 1: 486 raise RelaxError("Only one interatomic data container, hence dipole-dipole interaction, is supported per spin.") 487 488 # Spin title. 489 file.write("\nspin " + spin_id + "\n") 490 491 file.write('%-14s' % "constants") 492 file.write('%-6i' % res_num) 493 file.write('%-7s' % spin.isotope) 494 file.write('%-8.4f' % (periodic_table.gyromagnetic_ratio(spin.isotope) / 1e7)) 495 file.write('%-8.3f' % (interatoms[0].r * 1e10)) 496 file.write('%-8.3f\n' % (spin.csa * 1e6)) 497 498 file.write('%-10s' % "vector") 499 file.write('%-4s' % atom1) 500 file.write('%-4s\n' % atom2)

501 502

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

504 """Create the script 'run.sh' for the execution of Modelfree4. 505 506 @param file: The writable file object. 507 @type file: file object 508 @keyword binary: The name of the Modelfree4 binary file. This can include the path to the 509 binary. 510 @type binary: str 511 @keyword dir: The directory to copy the PDB file to. 512 @type dir: str 513 """ 514 515 # Check for the diffusion tensor. 516 if not hasattr(cdp, 'diff_tensor'): 517 raise RelaxNoTensorError('diffusion') 518 519 file.write("#! /bin/sh\n") 520 file.write(binary + " -i mfin -d mfdata -p mfpar -m mfmodel -o mfout -e out") 521 if cdp.diff_tensor.type != 'sphere': 522 # Copy the first pdb file to the model directory so there are no problems with the existance of *.rotate files. 523 mol = cdp.structure.structural_data[0].mol[0] 524 system('cp ' + mol.file_path + sep + mol.file_name + ' ' + dir) 525 file.write(" -s " + mol.file_name) 526 file.write("\n")

527 528

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

530 """Execute Modelfree4. 531 532 BUG: Control-C during execution causes the cwd to stay as dir. 533 534 535 @param dir: The optional directory where the script is located. 536 @type dir: str or None 537 @param force: A flag which if True will cause any pre-existing files to be overwritten by 538 Modelfree4. 539 @type force: bool 540 @param binary: The name of the Modelfree4 binary file. This can include the path to the 541 binary. 542 @type binary: str 543 """ 544 545 # Check for the diffusion tensor. 546 if not hasattr(cdp, 'diff_tensor'): 547 raise RelaxNoTensorError('diffusion') 548 549 # The current directory. 550 orig_dir = getcwd() 551 552 # The directory. 553 if dir == None: 554 dir = pipes.cdp_name() 555 if not access(dir, F_OK): 556 raise RelaxDirError('Modelfree4', dir) 557 558 # Change to this directory. 559 chdir(dir) 560 561 # Catch failures and return to the correct directory. 562 try: 563 # Python 2.3 and earlier. 564 if Popen == None: 565 raise RelaxError("The subprocess module is not available in this version of Python.") 566 567 # Test if the 'mfin' input file exists. 568 if not access('mfin', F_OK): 569 raise RelaxFileError('mfin input', 'mfin') 570 571 # Test if the 'mfdata' input file exists. 572 if not access('mfdata', F_OK): 573 raise RelaxFileError('mfdata input', 'mfdata') 574 575 # Test if the 'mfmodel' input file exists. 576 if not access('mfmodel', F_OK): 577 raise RelaxFileError('mfmodel input', 'mfmodel') 578 579 # Test if the 'mfpar' input file exists. 580 if not access('mfpar', F_OK): 581 raise RelaxFileError('mfpar input', 'mfpar') 582 583 # Test if the 'PDB' input file exists. 584 if cdp.diff_tensor.type != 'sphere': 585 pdb = cdp.structure.structural_data[0].mol[0].file_name 586 if not access(pdb, F_OK): 587 raise RelaxFileError('PDB', pdb) 588 else: 589 pdb = None 590 591 # Remove the file 'mfout' and '*.out' if the force flag is set. 592 if force: 593 for file in listdir(getcwd()): 594 if search('out$', file) or search('rotate$', file): 595 remove(file) 596 597 # Test the binary file string corresponds to a valid executable. 598 test_binary(binary) 599 600 # Execute Modelfree4. 601 if pdb: 602 cmd = binary + ' -i mfin -d mfdata -p mfpar -m mfmodel -o mfout -e out -s ' + pdb 603 else: 604 cmd = binary + ' -i mfin -d mfdata -p mfpar -m mfmodel -o mfout -e out' 605 pipe = Popen(cmd, shell=True, stdin=PIPE, stdout=PIPE, stderr=PIPE, close_fds=False) 606 607 # Close the pipe. 608 pipe.stdin.close() 609 610 # Write to stdout. 611 for line in pipe.stdout.readlines(): 612 # Decode Python 3 byte arrays. 613 if hasattr(line, 'decode'): 614 line = line.decode() 615 616 # Write. 617 sys.stdout.write(line) 618 619 # Write to stderr. 620 for line in pipe.stderr.readlines(): 621 # Decode Python 3 byte arrays. 622 if hasattr(line, 'decode'): 623 line = line.decode() 624 625 # Write. 626 sys.stderr.write(line) 627 628 # Failure. 629 except: 630 # Change back to the original directory. 631 chdir(orig_dir) 632 633 # Reraise the error. 634 raise 635 636 # Change back to the original directory. 637 chdir(orig_dir)

638 639

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

641 """Extract the Modelfree4 results out of the 'mfout' file. 642 643 @param dir: The directory containing the 'mfout' file. 644 @type dir: str or None 645 @keyword spin_id: The spin identification string. 646 @type spin_id: str or None 647 """ 648 649 # Test if sequence data is loaded. 650 if not exists_mol_res_spin_data(): 651 raise RelaxNoSequenceError 652 653 # Check for the diffusion tensor. 654 if not hasattr(cdp, 'diff_tensor'): 655 raise RelaxNoTensorError('diffusion') 656 657 # The directory. 658 if dir == None: 659 dir = pipes.cdp_name() 660 if not access(dir, F_OK): 661 raise RelaxDirError('Modelfree4', dir) 662 663 # Test if the file exists. 664 if not access(dir + sep+'mfout', F_OK): 665 raise RelaxFileError('Modelfree4', dir + sep+'mfout') 666 667 # Determine the parameter set. 668 model_type = determine_model_type() 669 670 # Open the file. 671 mfout_file = open(dir + sep+'mfout', 'r') 672 mfout_lines = mfout_file.readlines() 673 mfout_file.close() 674 675 # Get the section line positions of the mfout file. 676 global_chi2_pos, diff_pos, s2_pos, s2f_pos, s2s_pos, te_pos, rex_pos, chi2_pos = line_positions(mfout_lines) 677 678 # Find out if simulations were carried out. 679 sims = 0 680 for i in range(len(mfout_lines)): 681 if search('_iterations', mfout_lines[i]): 682 row = mfout_lines[i].split() 683 sims = int(row[1]) 684 685 # Global data. 686 if model_type in ['all', 'diff']: 687 # Global chi-squared. 688 row = mfout_lines[global_chi2_pos].split() 689 cdp.chi2 = float(row[1]) 690 691 # Spherical diffusion tensor. 692 if cdp.diff_tensor.type == 'sphere': 693 # Split the lines. 694 tm_row = mfout_lines[diff_pos].split() 695 696 # Set the params. 697 cdp.diff_tensor.set(param='tm', value=float(tm_row[2])) 698 699 # Spheroid diffusion tensor. 700 else: 701 # Split the lines. 702 tm_row = mfout_lines[diff_pos].split() 703 dratio_row = mfout_lines[diff_pos+1].split() 704 theta_row = mfout_lines[diff_pos+2].split() 705 phi_row = mfout_lines[diff_pos+3].split() 706 707 # Set the params. 708 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']) 709 710 # Loop over the sequence. 711 pos = 0 712 for spin, mol_name, res_num, res_name in spin_loop(spin_id, full_info=True): 713 # Skip deselected residues. 714 if not spin.select: 715 continue 716 717 # Get the residue number from the mfout file. 718 mfout_res_num = int(mfout_lines[s2_pos + pos].split()[0]) 719 720 # Skip the spin if the residue doesn't match. 721 if mfout_res_num != res_num: 722 continue 723 724 # Test that the model has been set (needed to differentiate between te and ts). 725 if not hasattr(spin, 'model'): 726 raise RelaxNoModelError 727 728 # Get the S2 data. 729 if 's2' in spin.params: 730 spin.s2, spin.s2_err = get_mf_data(mfout_lines, s2_pos + pos) 731 732 # Get the S2f data. 733 if 's2f' in spin.params or 's2s' in spin.params: 734 spin.s2f, spin.s2f_err = get_mf_data(mfout_lines, s2f_pos + pos) 735 736 # Get the S2s data. 737 if 's2f' in spin.params or 's2s' in spin.params: 738 spin.s2s, spin.s2s_err = get_mf_data(mfout_lines, s2s_pos + pos) 739 740 # Get the te data. 741 if 'te' in spin.params: 742 spin.te, spin.te_err = get_mf_data(mfout_lines, te_pos + pos) 743 spin.te = spin.te / 1e12 744 spin.te_err = spin.te_err / 1e12 745 746 # Get the ts data. 747 if 'ts' in spin.params: 748 spin.ts, spin.ts_err = get_mf_data(mfout_lines, te_pos + pos) 749 spin.ts = spin.ts / 1e12 750 spin.ts_err = spin.ts_err / 1e12 751 752 # Get the Rex data. 753 if 'rex' in spin.params: 754 spin.rex, spin.rex_err = get_mf_data(mfout_lines, rex_pos + pos) 755 spin.rex = spin.rex / (2.0 * pi * cdp.spectrometer_frq[cdp.ri_ids[0]])**2 756 spin.rex_err = spin.rex_err / (2.0 * pi * cdp.spectrometer_frq[cdp.ri_ids[0]])**2 757 758 # Get the chi-squared data. 759 if not sims: 760 row = mfout_lines[chi2_pos + pos].split() 761 spin.chi2 = float(row[1]) 762 else: 763 # 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). 764 row = mfout_lines[chi2_pos + 2 + 22*pos].split() 765 spin.chi2 = float(row[1]) 766 767 # Increment the residue position. 768 pos = pos + 1

769 770

771 -def get_mf_data(mfout_lines, pos):

772 """Extract the model-free data from the given position of the mfout file. 773 774 This method is designed to catch a number of bugs in Modelfree4's mfout file. 775 776 The first bug is the presence of a series of '*' characters causing a fusion of two columns. 777 This is handled by splitting by the '*' char and then returning the first element. 778 779 The second bug is when the floating point number is too big to fit into Modelfree4's string 780 format limit of 15.3f. This results in a results line such as: 781 782 246 10000.00019682363392.000 1 0.000 0.000 0.000 0.000 783 784 This is caught by scanning for two '.' characters in the column, and handled by assuming 785 that every floating point number will have three decimal characters. 786 787 @param mfout_lines: A list of all the lines of the mfout file. 788 @type mfout_lines: list of str 789 @param pos: The mfout line position. 790 @type pos: int 791 @return: The value and error. 792 @rtype: tuple of 2 floats 793 """ 794 795 # Split the line up. 796 row = mfout_lines[pos].split() 797 798 # The value and error, assuming a bug free mfout file. 799 val = row[1] 800 err = row[4] 801 802 # The Modelfree4 '*' column fusion bug. 803 if search('\*', val) or search('\*', err): 804 # Split by the '*' character. 805 val_row = val.split('*') 806 err_row = err.split('*') 807 808 # The value and error (the first elements). 809 val = val_row[0] 810 err = err_row[0] 811 812 # The Modelfree4 large float column fusion bug. 813 new_row = [] 814 fused = False 815 for element in row: 816 # Count the number of '.' characters. 817 num = element.count('.') 818 819 # Catch two or more '.' characters. 820 if num > 1: 821 # Set the fused flag. 822 fused = True 823 824 # Loop over each fused number. 825 for i in range(num): 826 # Find the index of the first '.'. 827 index = element.find('.') 828 829 # The first number (index + decimal point + 3 decimal chars). 830 new_row.append(element[0:index+4]) 831 832 # Strip the first number from the element for the next loop iteration. 833 element = element[index+4:] 834 835 # Otherwise the column element is fine. 836 else: 837 new_row.append(element) 838 839 # Bug has been caught. 840 if fused: 841 val = new_row[1] 842 err = new_row[4] 843 844 # Return the value and error, as floats. 845 return float(val), float(err)

846 847

848 -def line_positions(mfout_lines):

849 """Function for getting the section positions (line number) of the mfout file. 850 851 @param mfout_lines: A list of all the lines of the mfout file. 852 @type mfout_lines: list of str 853 @return: The line indices where the s2, s2f, s2s, te, rex, and chi2 sections 854 start in the mfout file. 855 @rtype: tuple of int 856 """ 857 858 # Loop over the file. 859 i = 0 860 while i < len(mfout_lines): 861 # Global chi2. 862 if match('data_chi_square', mfout_lines[i]): 863 global_chi2_pos = i + 1 864 865 # Diffusion tensor. 866 if match('data_diffusion_tensor', mfout_lines[i]): 867 diff_pos = i + 3 868 869 # Model-free data. 870 if match('data_model_1', mfout_lines[i]): 871 # Shift down two lines (to avoid the lines not starting with a space). 872 i = i + 2 873 874 # Walk through all the data. 875 while True: 876 # Break once the end of the data section is reached. 877 if not mfout_lines[i] == '\n' and not search('^ ', mfout_lines[i]): 878 break 879 880 # Split the line up. 881 row = mfout_lines[i].split() 882 883 # S2 position (skip the heading and move to the first residue). 884 if len(row) == 2 and row[0] == 'S2': 885 s2_pos = i + 1 886 887 # S2f position (skip the heading and move to the first residue). 888 if len(row) == 2 and row[0] == 'S2f': 889 s2f_pos = i + 1 890 891 # S2s position (skip the heading and move to the first residue). 892 if len(row) == 2 and row[0] == 'S2s': 893 s2s_pos = i + 1 894 895 # te position (skip the heading and move to the first residue). 896 if len(row) == 2 and row[0] == 'te': 897 te_pos = i + 1 898 899 # Rex position (skip the heading and move to the first residue). 900 if len(row) == 2 and row[0] == 'Rex': 901 rex_pos = i + 1 902 903 # Move to the next line number. 904 i = i + 1 905 906 # Chi-squared values. 907 if match('data_sse', mfout_lines[i]): 908 chi2_pos = i + 3 909 910 # Move to the next line number. 911 i = i + 1 912 913 # Return the positions. 914 return global_chi2_pos, diff_pos, s2_pos, s2f_pos, s2s_pos, te_pos, rex_pos, chi2_pos

915

Source Code for Module pipe_control.palmer