pipe_control.palmer

1 ############################################################################### 2 # # 3 # Copyright (C) 2001-2009,2011-2012,2017 Edward d'Auvergne # 4 # Copyright (C) 2008 Sebastien Morin # 5 # # 6 # This file is part of the program relax (http://www.nmr-relax.com). # 7 # # 8 # This program is free software: you can redistribute it and/or modify # 9 # it under the terms of the GNU General Public License as published by # 10 # the Free Software Foundation, either version 3 of the License, or # 11 # (at your option) any later version. # 12 # # 13 # This program is distributed in the hope that it will be useful, # 14 # but WITHOUT ANY WARRANTY; without even the implied warranty of # 15 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # 16 # GNU General Public License for more details. # 17 # # 18 # You should have received a copy of the GNU General Public License # 19 # along with this program. If not, see <http://www.gnu.org/licenses/>. # 20 # # 21 ############################################################################### 22 23 # Module docstring. 24 """Module for interfacing with Art Palmer's Modelfree 4 program.""" 25 26 27 # Dependencies. 28 import dep_check 29 30 # Python module imports. 31 from math import pi 32 from os import F_OK, access, chdir, chmod, getcwd, listdir, remove, sep, system 33 from re import match, search 34 import signal 35 from stat import S_IRWXU, S_IRGRP, S_IROTH 36 PIPE, Popen = None, None 37 if dep_check.subprocess_module: 38 from subprocess import PIPE, Popen 39 import sys 40 41 # relax module imports. 42 from lib.errors import RelaxError, RelaxDirError, RelaxFileError, RelaxNoInteratomError, RelaxNoModelError, RelaxNoPdbError, RelaxNoSequenceError, RelaxNoTensorError 43 from lib.io import mkdir_nofail, open_write_file, test_binary 44 from lib.periodic_table import periodic_table 45 from pipe_control import diffusion_tensor, pipes 46 from pipe_control.interatomic import return_interatom_list 47 from pipe_control.mol_res_spin import exists_mol_res_spin_data, spin_loop 48 from pipe_control.pipes import check_pipe 49 from specific_analyses.model_free.model import determine_model_type 50 51

52 -def __deselect_spins():

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

77 78

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

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

179 180

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

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

251 252

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

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

364 365

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

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

464 465

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

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

503 504

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

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

529 530

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

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

647 648

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

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

778 779

780 -def get_mf_data(mfout_lines, pos):

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

855 856

857 -def line_positions(mfout_lines):

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

924

Source Code for Module pipe_control.palmer