1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 """Module for interfacing with Dasha."""
24
25
26 import dep_check
27
28
29 from math import pi
30 from os import F_OK, access, chdir, getcwd, sep
31 PIPE, Popen = None, None
32 if dep_check.subprocess_module:
33 from subprocess import PIPE, Popen
34 import sys
35
36
37 from lib.errors import RelaxDirError, RelaxError, RelaxFileError, RelaxNoPdbError, RelaxNoSequenceError, RelaxNoTensorError
38 from lib.io import extract_data, mkdir_nofail, open_write_file, strip, test_binary
39 from pipe_control import angles, diffusion_tensor, pipes, value
40 from pipe_control.interatomic import return_interatom_list
41 from pipe_control.mol_res_spin import exists_mol_res_spin_data, first_residue_num, last_residue_num, residue_loop, return_spin, spin_loop
42 from pipe_control.spectrometer import loop_frequencies
43 from specific_analyses.setup import model_free_obj
44
45
71
72
73 -def create(algor='LM', dir=None, force=False):
126
127
129 """Create the Dasha script file.
130
131 @param file: The opened file descriptor.
132 @type file: file object
133 @param model_type: The model-free model type.
134 @type model_type: str
135 @param algor: The optimisation algorithm to use. This can be the Levenberg-Marquardt algorithm 'LM' or the Newton-Raphson algorithm 'NR'.
136 @type algor: str
137 """
138
139
140 file.write("# Delete all data.\n")
141 file.write("del 1 10000\n")
142
143
144 file.write("\n# Nucleus type.\n")
145 nucleus = None
146 for spin in spin_loop():
147
148 if spin.isotope == '1H':
149 continue
150
151
152 if nucleus and spin.isotope != nucleus:
153 raise RelaxError("The nuclei '%s' and '%s' do not match, relax can only handle one nucleus type in Dasha." % (nucleus, spin.isotope))
154
155
156 if not nucleus:
157 nucleus = spin.isotope
158
159
160 if nucleus == '15N':
161 nucleus = 'N15'
162 elif nucleus == '13C':
163 nucleus = 'C13'
164 else:
165 raise RelaxError("Cannot handle the nucleus type '%s' within Dasha." % nucleus)
166 file.write("set nucl %s\n" % nucleus)
167
168
169 file.write("\n# Number of frequencies.\n")
170 file.write("set n_freq %s\n" % cdp.spectrometer_frq_count)
171
172
173 file.write("\n# Frequency values.\n")
174 count = 1
175 for frq in loop_frequencies():
176 file.write("set H1_freq %s %s\n" % (frq / 1e6, count))
177 count += 1
178
179
180 file.write("\n# Set the diffusion tensor.\n")
181 if model_type != 'local_tm':
182
183 if cdp.diff_tensor.type == 'sphere':
184 file.write("set tr %s\n" % (cdp.diff_tensor.tm / 1e-9))
185
186
187 elif cdp.diff_tensor.type == 'spheroid':
188 file.write('set tr %s\n' % (cdp.diff_tensor.tm / 1e-9))
189
190
191 elif cdp.diff_tensor.type == 'ellipsoid':
192
193 Dx, Dy, Dz = diffusion_tensor.return_eigenvalues()
194
195
196 file.write("set tr %s\n" % (cdp.diff_tensor.tm / 1e-9))
197 file.write("set D1/D3 %s\n" % (Dx / Dz))
198 file.write("set D2/D3 %s\n" % (Dy / Dz))
199
200
201 file.write("set alfa %s\n" % (cdp.diff_tensor.alpha / (2.0 * pi) * 360.0))
202 file.write("set betta %s\n" % (cdp.diff_tensor.beta / (2.0 * pi) * 360.0))
203 file.write("set gamma %s\n" % (cdp.diff_tensor.gamma / (2.0 * pi) * 360.0))
204
205
206 file.write("\n# Reading the relaxation data.\n")
207 file.write("echo Reading the relaxation data.\n")
208 noe_index = 1
209 r1_index = 1
210 r2_index = 1
211 for ri_id in cdp.ri_ids:
212
213 if cdp.ri_type[ri_id] == 'NOE':
214
215 number = noe_index
216
217
218 data_type = 'noe'
219
220
221 noe_index = noe_index + 1
222
223
224 elif cdp.ri_type[ri_id] == 'R1':
225
226 number = r1_index
227
228
229 data_type = '1/T1'
230
231
232 r1_index = r1_index + 1
233
234
235 elif cdp.ri_type[ri_id] == 'R2':
236
237 number = r2_index
238
239
240 data_type = '1/T2'
241
242
243 r2_index = r2_index + 1
244
245
246 if number == 1:
247 file.write("\nread < %s\n" % data_type)
248 else:
249 file.write("\nread < %s %s\n" % (data_type, number))
250
251
252 for residue in residue_loop():
253
254 spin = residue.spin[0]
255
256
257 if not spin.select:
258 continue
259
260
261 if len(spin.ri_data) != len(cdp.ri_ids) or spin.ri_data[ri_id] == None:
262 spin.select = False
263 continue
264
265
266 file.write("%s %s %s\n" % (residue.num, spin.ri_data[ri_id], spin.ri_data_err[ri_id]))
267
268
269 file.write("exit\n")
270
271
272 if model_type == 'mf':
273
274 for residue in residue_loop():
275
276 spin = residue.spin[0]
277
278
279 if not spin.select:
280 continue
281
282
283 interatoms = return_interatom_list(spin._spin_ids[0])
284 if len(interatoms) == 0:
285 raise RelaxNoInteratomError
286 elif len(interatoms) > 1:
287 raise RelaxError("Only one interatomic data container, hence dipole-dipole interaction, is supported per spin.")
288
289
290 file.write("\n\n\n# Residue %s\n\n" % residue.num)
291
292
293 file.write("echo Optimisation of residue %s\n" % residue.num)
294
295
296 file.write("\n# Select the residue.\n")
297 file.write("set cres %s\n" % residue.num)
298
299
300 if cdp.diff_tensor.type == 'spheroid':
301 file.write("set teta %s\n" % spin.alpha)
302
303
304 elif cdp.diff_tensor.type == 'ellipsoid':
305 file.write("\n# Setting the spherical angles of the XH vector in the ellipsoid diffusion frame.\n")
306 file.write("set teta %s\n" % spin.theta)
307 file.write("set fi %s\n" % spin.phi)
308
309
310 if 'ts' in spin.params:
311 jmode = 3
312 elif 'te' in spin.params:
313 jmode = 2
314 elif 's2' in spin.params:
315 jmode = 1
316
317
318 if 'rex' in spin.params:
319 exch = True
320 else:
321 exch = False
322
323
324 if cdp.diff_tensor.type == 'sphere':
325 anis = False
326 else:
327 anis = True
328
329
330 if cdp.diff_tensor.type == 'spheroid':
331 sym = True
332 else:
333 sym = False
334
335
336 file.write("\n# Set the jmode.\n")
337 file.write("set def jmode %s" % jmode)
338 if exch:
339 file.write(" exch")
340 if anis:
341 file.write(" anis")
342 if sym:
343 file.write(" sym")
344 file.write("\n")
345
346
347 file.write("\n# Parameter default values.\n")
348 file.write("reset jmode %s\n" % residue.num)
349
350
351 file.write("\n# Bond length.\n")
352 file.write("set r_hx %s\n" % (interatoms[0].r / 1e-10))
353
354
355 file.write("\n# CSA value.\n")
356 file.write("set csa %s\n" % (spin.csa / 1e-6))
357
358
359 if not 'tf' in spin.params and jmode == 3:
360 file.write("\n# Fix the tf parameter.\n")
361 file.write("fix tf 0\n")
362
363
364 file.write("\n\n\n# Optimisation of all residues.\n")
365 if algor == 'LM':
366 file.write("lmin %s %s" % (first_residue_num(), last_residue_num()))
367 elif algor == 'NR':
368 file.write("min %s %s" % (first_residue_num(), last_residue_num()))
369
370
371 file.write("\n# Show the results.\n")
372 file.write("echo\n")
373 file.write("show all\n")
374
375
376 file.write("\n# Write the results.\n")
377 file.write("write s2.out S\n")
378 file.write("write s2f.out Sf\n")
379 file.write("write s2s.out Ss\n")
380 file.write("write te.out te\n")
381 file.write("write tf.out tf\n")
382 file.write("write ts.out ts\n")
383 file.write("write rex.out rex\n")
384 file.write("write chi2.out F\n")
385
386 else:
387 raise RelaxError("Optimisation of the parameter set '%s' currently not supported." % model_type)
388
389
391 """Execute Dasha.
392
393 @param dir: The optional directory where the script is located.
394 @type dir: str or None
395 @param force: A flag which if True will cause any pre-existing files to be overwritten by Dasha.
396 @type force: bool
397 @param binary: The name of the Dasha binary file. This can include the path to the binary.
398 @type binary: str
399 """
400
401
402 test_binary(binary)
403
404
405 orig_dir = getcwd()
406
407
408 if dir == None:
409 dir = pipes.cdp_name()
410 if not access(dir, F_OK):
411 raise RelaxDirError('Dasha', dir)
412
413
414 chdir(dir)
415
416
417 try:
418
419 if not access('dasha_script', F_OK):
420 raise RelaxFileError('dasha script', 'dasha_script')
421
422
423 if Popen == None:
424 raise RelaxError("The subprocess module is not available in this version of Python.")
425
426
427 pipe = Popen(binary, shell=True, stdin=PIPE, stdout=PIPE, stderr=PIPE, close_fds=False)
428
429
430 script = open('dasha_script')
431 lines = script.readlines()
432 script.close()
433 for line in lines:
434
435 if hasattr(line, 'encode'):
436 line = line.encode()
437
438
439 pipe.stdin.write(line)
440
441
442 pipe.stdin.close()
443
444
445 for line in pipe.stdout.readlines():
446
447 if hasattr(line, 'decode'):
448 line = line.decode()
449
450
451 sys.stdout.write(line)
452
453
454 for line in pipe.stderr.readlines():
455
456 if hasattr(line, 'decode'):
457 line = line.decode()
458
459
460 sys.stderr.write(line)
461
462
463 except:
464
465 chdir(orig_dir)
466
467
468 raise
469
470
471 chdir(orig_dir)
472
473
474 sys.stdout.write("\n\n")
475
476
478 """Extract the data from the Dasha results files.
479
480 @param dir: The optional directory where the results file is located.
481 @type dir: str or None
482 """
483
484
485 if not exists_mol_res_spin_data():
486 raise RelaxNoSequenceError
487
488
489 if dir == None:
490 dir = pipes.cdp_name()
491 if not access(dir, F_OK):
492 raise RelaxDirError('Dasha', dir)
493
494
495 for param in ['s2', 's2f', 's2s', 'te', 'tf', 'ts', 'rex']:
496
497 file_name = dir + sep + param + '.out'
498
499
500 if not access(file_name, F_OK):
501 raise RelaxFileError('Dasha', file_name)
502
503
504 if param in ['te', 'tf', 'ts']:
505 scaling = 1e-9
506 elif param == 'rex':
507 scaling = 1.0 / (2.0 * pi * cdp.spectrometer_frq[cdp.ri_ids[0]]) ** 2
508 else:
509 scaling = 1.0
510
511
512 data = read_results(file=file_name, scaling=scaling)
513
514
515 for i in range(len(data)):
516 value.set(val=data[i][1], param=param, spin_id=data[i][0])
517 value.set(val=data[i][0], param=param, spin_id=data[i][0], error=True)
518
519
520 for spin in spin_loop():
521
522 if param in spin.params:
523 continue
524
525
526 setattr(spin, param.lower(), None)
527
528
529 file_name = dir + sep+'chi2.out'
530
531
532 if not access(file_name, F_OK):
533 raise RelaxFileError('Dasha', file_name)
534
535
536 data = read_results(file=file_name)
537
538
539 for i in range(len(data)):
540 spin = return_spin(data[i][0])
541 spin.chi2 = data[i][1]
542
543
545 """Extract the data from the Dasha results file.
546
547 @keyword file: The name of the file to open.
548 @type file: str
549 @keyword dir: The directory containing the file (defaults to the current directory if None).
550 @type dir: str or None
551 @keyword scaling: The parameter scaling factor.
552 @type scaling: float
553 """
554
555
556 data = extract_data(file=file, dir=dir)
557
558
559 data = strip(data)
560
561
562 new_data = []
563 for i in range(len(data)):
564 spin_id = ':%s@N' % data[i][0]
565 value = float(data[i][1]) * scaling
566 error = float(data[i][2]) * scaling
567 new_data.append([spin_id, value, error])
568
569
570 return new_data
571