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.physical_constants import return_gyromagnetic_ratio 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 specific_analyses.model_free.model import determine_model_type 47 48

49 -def __deselect_spins():

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

74 75

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

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

176 177

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

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

248 249

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

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

361 362

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

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

461 462

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

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

500 501

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

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

526 527

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

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

637 638

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

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

768 769

770 -def get_mf_data(mfout_lines, pos):

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

845 846

847 -def line_positions(mfout_lines):

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

914

Source Code for Module pipe_control.palmer