specific_analyses.relax

1 ############################################################################### 2 # # 3 # Copyright (C) 2013-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 """Functions for interfacing with Flemming Hansen's CATIA program.""" 24 25 # Dependencies. 26 import dep_check 27 28 # Python module imports. 29 from os import F_OK, access, chdir, getcwd, sep 30 PIPE, Popen = None, None 31 if dep_check.subprocess_module: 32 from subprocess import PIPE, Popen 33 import sys 34 35 # relax module imports. 36 from lib.errors import RelaxError, RelaxDirError 37 from lib.io import mkdir_nofail, open_write_file, test_binary 38 from pipe_control import pipes 39 from pipe_control.mol_res_spin import check_mol_res_spin_data, spin_loop 40 from specific_analyses.relax_disp.checks import check_model_type, check_spectra_id_setup 41 from specific_analyses.relax_disp.data import loop_frq, loop_offset_point, return_param_key_from_data 42 from specific_analyses.relax_disp.variables import EXP_TYPE_CPMG_SQ 43 44

45 -def catia_execute(file='Fit.catia', dir=None, binary=None):

46 """Create the CATIA input files. 47 48 @keyword file: The main CATIA execution file. 49 @type file: str 50 @keyword dir: The optional directory to place the files into. If None, then the files will be placed into the current directory. 51 @type dir: str or None 52 @keyword binary: The name of the CATIA binary file. This can include the path to the binary. 53 @type binary: str 54 """ 55 56 # The current directory. 57 orig_dir = getcwd() 58 59 # Test the binary file string corresponds to a valid executable. 60 test_binary(binary) 61 62 # The directory. 63 if dir == None: 64 dir = orig_dir 65 if not access(dir, F_OK): 66 raise RelaxDirError('CATIA', dir) 67 68 # Change to the directory with the CATIA input files. 69 chdir(dir) 70 71 # Catch failures and return to the correct directory. 72 try: 73 # Execute CATIA. 74 cmd = binary + ' < Fit.catia' 75 pipe = Popen(cmd, shell=True, stdin=PIPE, stdout=PIPE, stderr=PIPE, close_fds=False) 76 77 # Close the pipe. 78 pipe.stdin.close() 79 80 # Write to stdout. 81 for line in pipe.stdout.readlines(): 82 # Decode Python 3 byte arrays. 83 if hasattr(line, 'decode'): 84 line = line.decode() 85 86 # Write. 87 sys.stdout.write(line) 88 89 # Write to stderr. 90 for line in pipe.stderr.readlines(): 91 # Decode Python 3 byte arrays. 92 if hasattr(line, 'decode'): 93 line = line.decode() 94 95 # Write. 96 sys.stderr.write(line) 97 98 # Failure. 99 except: 100 # Change back to the original directory. 101 chdir(orig_dir) 102 103 # Reraise the error. 104 raise 105 106 # Change back to the original directory. 107 chdir(orig_dir)

108 109

110 -def catia_input(file='Fit.catia', dir=None, output_dir='output', force=False):

111 """Create the CATIA input files. 112 113 @keyword file: The main CATIA execution file. 114 @type file: str 115 @keyword dir: The optional directory to place the files into. If None, then the files will be placed into the current directory. 116 @type dir: str or None 117 @keyword output_dir: The CATIA output directory, located within the directory specified by the dir argument. This directory will be created. 118 @type output_dir: str 119 @keyword force: A flag which if True will cause all pre-existing files to be overwritten. 120 @type force: Bool 121 """ 122 123 # Data checks. 124 pipes.test() 125 check_mol_res_spin_data() 126 check_spectra_id_setup() 127 check_model_type() 128 129 # Check that this is CPMG data. 130 for id in cdp.spectrum_ids: 131 if cdp.exp_type[id] != 'SQ CPMG': 132 raise RelaxError("Only CPMG type data is supported.") 133 134 # Directory creation. 135 if dir != None: 136 mkdir_nofail(dir, verbosity=0) 137 138 # Create the R2eff files. 139 write_r2eff_files(input_dir='input_r2eff', base_dir=dir, force=force) 140 141 # Create the parameter files. 142 write_param_files(global_file="ParamGlobal.inp", set_file="ParamSet1.inp", dir=dir, force=force) 143 144 # Create the main execution file. 145 write_main_file(file=file, dir=dir, output_dir=output_dir, force=force) 146 147 # Create the output directory as needed by CATIA (it does not create it itself). 148 mkdir_nofail(dir + sep + output_dir, verbosity=0)

149 150

151 -def write_main_file(file=None, dir=None, output_dir=None, f_tol=1e-25, max_iter=10000000, r1=False, force=False):

152 """Create the main CATIA execution file. 153 154 @keyword file: The main CATIA execution file. 155 @type file: str 156 @keyword dir: The directory to place the files into. 157 @type dir: str or None 158 @keyword output_dir: The CATIA output directory, located within the directory specified by the dir argument. This directory will be created. 159 @type output_dir: str 160 @keyword r1: A flag which if True will cause the R1 data to be used for off-resonance effects. 161 @type r1: bool 162 @keyword force: A flag which if True will cause a pre-existing file to be overwritten. 163 @type force: bool 164 """ 165 166 # The file. 167 catia_in = open_write_file(file_name=file, dir=dir, force=force) 168 169 # The R2eff input sets. 170 for frq in loop_frq(): 171 frq_label = int(frq*1e-6) 172 file_name = "data_set_%i.inp" % frq_label 173 catia_in.write("ReadDataset(%s)\n" % file_name) 174 175 # Write out the data. 176 catia_in.write("ReadParam_Global(ParamGlobal.inp)\n") 177 catia_in.write("ReadParam_Local(ParamSet1.inp)\n") 178 catia_in.write("\n") 179 180 # The R1 data for off-resonance effects. 181 if r1: 182 catia_in.write("ReadParam(Omega;%s;0;1)\n" % shift_file) 183 for frq in loop_frq(): 184 frq_label = int(frq*1e-6) 185 param = "R1iph_%s" % frq_label 186 r1_file = "R1_%s.out" % frq_label 187 catia_in.write("ReadParam(%s;%s;0;1)\n" % (param, r1_file)) 188 catia_in.write("\n") 189 190 # Fix these off-resonance parameters. 191 catia_in.write("FreeLocalParam(all;Omega;false)\n") 192 for frq in loop_frq(): 193 frq_label = int(frq*1e-6) 194 param = "R1iph_%s" % frq_label 195 catia_in.write("FreeLocalParam(all;%s;false)\n" % param) 196 catia_in.write("\n") 197 198 # Minimisation. 199 catia_in.write("Minimize(print=y;tol=%s;maxiter=%i)\n" % (f_tol, max_iter)) 200 catia_in.write("\n") 201 202 # Plotting. 203 catia_in.write("PrintParam(%s/GlobalParam.fit;global)\n" % output_dir) 204 catia_in.write("PrintParam(%s/DeltaOmega.fit;DeltaO)\n" % output_dir) 205 catia_in.write("PrintData(%s/)\n" % output_dir) 206 catia_in.write("\n") 207 208 # Calculate the chi-squared value (not sure why, it's calculated in the minimisation). 209 catia_in.write("ChiSq(all;all)\n") 210 211 # Exit the program. 212 catia_in.write("exit(0)\n") 213 214 # Close the file. 215 catia_in.close()

216 217

218 -def write_param_files(global_file=None, kex=1500.0, pA=0.95, set_file=None, dir=None, r1=False, force=False):

219 """Create the CATIA parameter files. 220 221 @keyword global_file: The name of the global parameter file. 222 @type global_file: str 223 @keyword set_file: The name of the parameter set file. 224 @type set_file: str 225 @keyword dir: The base directory to place the files into. 226 @type dir: str 227 @keyword r1: A flag which if True will cause the R1 data to be used for off-resonance effects. 228 @type r1: bool 229 @keyword force: A flag which if True will cause a pre-existing file to be overwritten. 230 @type force: bool 231 """ 232 233 # Open the global parameter file. 234 param_file = open_write_file(file_name=global_file, dir=dir, force=force) 235 236 # Set the starting values. 237 param_file.write("kex = %s\n" % kex) 238 param_file.write("pb = %s\n" % (1.0-pA)) 239 240 # Close the file. 241 param_file.close() 242 243 # Open the 1st parameter set file. 244 set_file = open_write_file(file_name=set_file, dir=dir, force=force) 245 246 # The parameter format. 247 params = ['Delta0'] 248 values = [0.5] 249 if r1: 250 for frq in loop_frq(): 251 params.append("R1iph_%s" % frq_label(frq)) 252 values.append(1.5) 253 for frq in loop_frq(): 254 params.append("R0iph_%s" % frq_label(frq)) 255 values.append(5.0) 256 if r1: 257 params.append("Omega") 258 values.append(0.0) 259 260 # Write out the format. 261 set_file.write("format = (") 262 for i in range(len(params)): 263 if i != 0: 264 set_file.write(';') 265 set_file.write(params[i]) 266 set_file.write(")\n") 267 268 # Um?!? The average values maybe? 269 set_file.write("* = (") 270 for i in range(len(values)): 271 if i != 0: 272 set_file.write(';') 273 set_file.write("%s" % values[i]) 274 set_file.write(")\n") 275 276 # Close the file. 277 set_file.close()

278 279

280 -def frq_label(frq):

281 """Convert the frequency in Hz to a label in MHz. 282 283 @param frq: The frequency to convert. 284 @type frq: float 285 @return: The frequency in MHz as a label. 286 @rtype: str 287 """ 288 289 # Convert and return. 290 return str(int(frq*1e-6))

291 292

293 -def write_r2eff_files(input_dir=None, base_dir=None, force=False):

294 """Create the CATIA R2eff input files. 295 296 @keyword input_dir: The special directory for the R2eff input files. 297 @type input_dir: str 298 @keyword base_dir: The base directory to place the files into. 299 @type base_dir: str 300 @keyword force: A flag which if True will cause a pre-existing file to be overwritten. 301 @type force: bool 302 """ 303 304 # Create the directory for the R2eff files for each field and spin. 305 dir = base_dir + sep + input_dir 306 mkdir_nofail(dir, verbosity=0) 307 308 # Determine the isotope information. 309 isotope = None 310 for spin in spin_loop(skip_desel=True): 311 if hasattr(spin, 'isotope'): 312 if isotope == None: 313 isotope = spin.isotope 314 elif spin.isotope != isotope: 315 raise RelaxError("CATIA only supports one spin type.") 316 if isotope == None: 317 raise RelaxError("The spin isotopes have not been specified.") 318 319 # Isotope translation. 320 if isotope == '1H': 321 isotope = 'H1' 322 elif isotope == '13C': 323 isotope = 'C13' 324 elif isotope == '15N': 325 isotope = 'N15' 326 327 # Loop over the frequencies. 328 for frq, mi in loop_frq(return_indices=True): 329 # The frequency string in MHz. 330 frq_string = int(frq*1e-6) 331 332 # The set files. 333 file_name = "data_set_%i.inp" % frq_string 334 set_file = open_write_file(file_name=file_name, dir=base_dir, force=force) 335 id = frq_string 336 set_file.write("ID=%s\n" % id) 337 set_file.write("Sfrq = %s\n" % frq_string) 338 set_file.write("Temperature = %s\n" % 0.0) 339 set_file.write("Nucleus = %s\n" % isotope) 340 set_file.write("Couplednucleus = %s\n" % 'H1') 341 set_file.write("Time_equil = %s\n" % 0.0) 342 set_file.write("Pwx_cp = %s\n" % 0.0) 343 set_file.write("Taub = %s\n" % 0.0) 344 set_file.write("Time_T2 = %s\n"% cdp.relax_time_list[0]) 345 set_file.write("Xcar = %s\n" % 0.0) 346 set_file.write("Seqfil = %s\n" % 'CW_CPMG') 347 set_file.write("Minerror = %s\n" % "(2.%;0.5/s)") 348 set_file.write("Basis = (%s)\n" % "Iph_7") 349 set_file.write("Format = (%i;%i;%i)\n" % (0, 1, 2)) 350 set_file.write("DataDirectory = %s\n" % (dir+sep)) 351 set_file.write("Data = (\n") 352 353 # Loop over the spins. 354 for spin, mol_name, res_num, res_name, spin_id in spin_loop(full_info=True, return_id=True, skip_desel=True): 355 # The file. 356 file_name = "spin%s_%i.cpmg" % (spin_id.replace('#', '_').replace(':', '_').replace('@', '_'), frq_string) 357 spin_file = open_write_file(file_name=file_name, dir=dir, force=force) 358 359 # Write the header. 360 spin_file.write("# %18s %20s %20s\n" % ("nu_cpmg(Hz)", "R2(1/s)", "Esd(R2)")) 361 362 # Loop over the dispersion points. 363 for offset, point, oi, di in loop_offset_point(exp_type=EXP_TYPE_CPMG_SQ, frq=frq, return_indices=True): 364 # The key. 365 key = return_param_key_from_data(exp_type=EXP_TYPE_CPMG_SQ, frq=frq, offset=offset, point=point) 366 367 # No data. 368 if key not in spin.r2eff: 369 continue 370 371 # Write out the data. 372 spin_file.write("%20.15f %20.15f %20.15f\n" % (point, spin.r2eff[key], spin.r2eff_err[key])) 373 374 # Close the file. 375 spin_file.close() 376 377 # Add the file name to the set. 378 catia_spin_id = "%i%s" % (res_num, spin.name) 379 set_file.write(" [%s;%s];\n" % (catia_spin_id, file_name)) 380 381 # Terminate the set file. 382 set_file.write(")\n") 383 set_file.close()

384

Source Code for Module specific_analyses.relax_disp.catia