File: prodrg_import.py

package info (click to toggle)
coot 1.1.18%2Bdfsg-3
  • links: PTS, VCS
  • area: main
  • in suites: sid
  • size: 219,964 kB
  • sloc: cpp: 495,934; python: 35,043; ansic: 26,143; lisp: 22,768; sh: 13,186; makefile: 2,746; awk: 441; xml: 245; csh: 14
file content (493 lines) | stat: -rw-r--r-- 20,028 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
cprodrg = "cprodrg"
# we cannot use full path as cprodrg is spawning refmac...
#cprodrg = "c:/Programs/CCP4-Packages/ccp4-6.1.13/bin/cprodrg.exe"

# we need $CLIBD to run prodrg (for prodrg.param), so check for it:
if not os.getenv("CLIBD"):
    # print "BL INFO:: potential problem to run prodrg (no CLIBD set)"
    # print "          try to fix this"
    bin_dir = os.path.dirname(cprodrg)
    base_dir = os.path.dirname(bin_dir)
    clibd = os.path.join(base_dir, "lib", "data")
    prodrg_params = os.path.join(clibd, "prodrg.param")
    if os.path.isfile(prodrg_params):
        os.environ["CLIBD"] = clibd
        print "BL INFO:: found prodrg.param"
        print "          and set CLIBD to", clibd
    else:
        # print "BL ERROR:: no prodrg.param available"
        # print "searched in ", clibd
        if coot.enhanced_ligand_coot_p():
            # we don't want to know that you can't find old programs
            pass
        else:
           print "Sorry, prodrg module not available"
        # should load rest?
        # i.e. stop here?
    

# these are the files that mdl-latest-time and sbase time-out functions
# look at.
# 
# if there is a prodrg-xyzin set the current-time to its mtime, else False
#
global prodrg_xyzin
global sbase_to_coot_tlc
prodrg_xyzin      = "prodrg-in.mdl"
sbase_to_coot_tlc = ".sbase-to-coot-comp-id"

# this is for BL win machine
home = os.getenv("HOME")
if (not home and coot_utils.is_windows()):
    home = os.getenv("COOT_HOME")
if home:
    prodrg_xyzin      = os.path.join(home, "Projects",
                                     "build-xp-python", "lbg", "prodrg-in.mdl")
else:
    print "BL DEBUG:: buggery home is", home  # FIXME
sbase_to_coot_tlc = "../../build-xp-python/lbg/.sbase-to-coot-comp-id"

# this is latest!!!!
prodrg_xyzin      = "prodrg-in.mdl"
sbase_to_coot_tlc = ".sbase-to-coot-comp-id"

# to be over-ridden by your favourite 3d conformer generator, if you like...
#
def import_from_3d_generator_from_mdl(file_name, comp_id):
    import_from_prodrg("mini-no", comp_id)


def import_ligand_with_overlay(prodrg_xyzout, prodrg_cif):
    
    # OK, so here we read the PRODRG files and
    # manipulate them.  We presume that the active
    # residue is quite like the input ligand from
    # prodrg.
    #
    # Read in the lib and coord output of PRODRG.  Then
    # overalay the new ligand onto the active residue
    # (just so that we can see it approximately
    # oriented). Then match the torsions from the new
    # ligand to the those of the active residue.  Then
    # overlay again so that we have the best match.
    #
    # We want to see just one molecule with the protein
    # and the new ligand.
    # add_ligand_delete_residue_copy_molecule provides
    # that for us.  We just colour it and undisplay the
    # other molecules.
    
    # overlap the imol_ligand residue if there are restraints for the
    # reference residue/ligand.
    #
    # Don't overlap if the reference residue/ligand is not a het-group.
    #
    # return overlapped status
    #
    def overlap_ligands_maybe(imol_ligand, imol_ref, chain_id_ref, res_no_ref):

        # we don't want to overlap-ligands if there is no dictionary
        # for the residue to be matched to.
        res_name = coot.residue_name(imol_ref, chain_id_ref, res_no_ref, "")
        restraints = monomer_restraints(res_name)
        if (not restraints):
            return False
        else:
            if not coot_utils.residue_has_hetatms_qm(imol_ref, chain_id_ref, res_no_ref, ""):
                return False
            else:
                print "----------- overlap-ligands %s %s %s %s ------------" \
                      %(imol_ligand, imol_ref, chain_id_ref, res_no_ref)
                # this can return the rtop operator or False (for fail of course).
                return overlap_ligands(imol_ligand, imol_ref, chain_id_ref, res_no_ref)
            
    # return the new molecule number.
    #
    def read_and_regularize(prodrg_xyzout):
        imol = coot.handle_read_draw_molecule_and_move_molecule_here(prodrg_xyzout)
        # speed up the minisation (and then restore setting).
        # No need to put the refinement step stuff in auto accept!?
        s = coot.dragged_refinement_steps_per_frame()
        coot.set_dragged_refinement_steps_per_frame(500)
        with AutoAccept():
            coot.regularize_residues( imol, [["", 1, ""]])
        coot.set_dragged_refinement_steps_per_frame(s)
        return imol

    # return the new molecule number
    # (only works with aa_ins_code of "")
    #
    def read_regularize_and_match_torsions_maybe(prodrg_xyzout, imol_ref,
                                                 chain_id_ref, res_no_ref):
        imol = coot.handle_read_draw_molecule_and_move_molecule_here(prodrg_xyzout)

        if (not have_restraints_for_qm(coot.residue_name(imol_ref, chain_id_ref,
                                                    res_no_ref, ""))):
            return False
        else:
            overlap_status = overlap_ligands_maybe(imol, imol_ref,
                                                   chain_id_ref, res_no_ref)
            
            # speed up the minisation (and then restore setting).
            # No need to put refinement steps in auto accept!?
            s = coot.dragged_refinement_steps_per_frame()
            coot.set_dragged_refinement_steps_per_frame(600)
            with AutoAccept():
                coot.regularize_residues(imol, [["", 1, ""]])
            coot.set_dragged_refinement_steps_per_frame(s)
            if overlap_status:
                coot.match_ligand_torsions(imol, imol_ref, chain_id_ref, res_no_ref)
        return imol

    # return True or False
    #
    def have_restraints_for_qm(res_name):
        restraints = monomer_restraints(res_name)
        return not restraints == False   # twisted but short

    # main line

    coot.read_cif_dictionary(prodrg_cif)

    # we do different things depending on whether or
    # not there is an active residue.  We need to test
    # for having an active residue here (currently we
    # presume that there is).
    # 
    # Similarly, if the aa_ins_code is non-null, let's
    # presume that we do not have an active residue.

    active_atom = active_residue()
    if ((not active_atom) or
        (not active_atom[3] == "")):  # aa_ins_code

        # then there is no active residue to match to
        read_and_regularize(prodrg_xyzout)
        # BL says:: maybe we should merge the ligand into
        # the protein molecule?! But what is the protein?
        # merge_molecules([active_atom[0]], imol)
    else:
        # we have an active residue to match to
        with coot_utils.UsingActiveAtom() as [aa_imol, aa_chain_id, aa_res_no, aa_ins_code,
                                   aa_atom_name, aa_alt_conf]:
            if not coot_utils.residue_is_close_to_screen_centre_qm(
                aa_imol, aa_chain_id, aa_res_no, ""):

                # not close, no overlap
                #
                read_and_regularize(prodrg_xyzout)
            else:
                # try overlap
                # BL says:: this overwrites 
                imol = read_regularize_and_match_torsions_maybe(prodrg_xyzout,
                                                                aa_imol, aa_chain_id, aa_res_no)

                overlapped_flag = overlap_ligands_maybe(imol, aa_imol,
                                                        aa_chain_id, aa_res_no)

                if overlapped_flag:
                    print "------ overlapped-flag was true!!!!!"
                    coot.set_mol_displayed(aa_imol, 0)
                    coot.set_mol_active(aa_imol, 0)
                    col = coot.get_molecule_bonds_colour_map_rotation(aa_imol)
                    new_col = col + 5 # a tiny amount
                    # new ligand specs, then "reference" ligand (to be deleted)
                    imol_replaced = coot.add_ligand_delete_residue_copy_molecule(
                        imol, "", 1,
                        aa_imol, aa_chain_id, aa_res_no)


                    coot.set_molecule_bonds_colour_map_rotation(imol_replaced, new_col)
                    coot.set_mol_displayed(imol, 0)
                    coot.set_mol_active(imol, 0)
                    coot.graphics_draw()

    return True
    # return False otherwise? When? FIXME

    
def import_from_prodrg(minimize_mode, res_name):

    import operator
    global prodrg_xyzin

    # main line of import_from_prodrg
    #
    prodrg_dir = "coot-ccp4"

    coot.make_directory_maybe(prodrg_dir)
    prodrg_xyzout = os.path.join(prodrg_dir, "prodrg-" + res_name + ".pdb")
    prodrg_cif    = os.path.join(prodrg_dir, "prodrg-out.cif")
    prodrg_log    = os.path.join(prodrg_dir, "prodrg.log")

    # FIXME? maybe, we could use 'enum' for mini-no rather than string?!
    # see news_status stuff
    # requires python >= 2.5 (shall we test?)
    mini_mode = "NO" if (minimize_mode == 'mini-no') else "PREP"
    # see if we have cprodrg
    if not (os.path.isfile(cprodrg) or
            coot_utils.command_in_path_qm(cprodrg)):
        coot.info_dialog("BL INFO:: No cprodrg found")
    else:
        status = coot_utils.popen_command(cprodrg,
                               ["XYZIN",  prodrg_xyzin,
                                "XYZOUT", prodrg_xyzout,
                                "LIBOUT", prodrg_cif],
                               ["MINI " + mini_mode, "END"],
                               prodrg_log, True)
        if isinstance(status, int):
            if (status == 0):
                import_ligand_with_overlay(prodrg_xyzout, prodrg_cif)


def get_file_latest_time(file_name):
    if not os.path.isfile(file_name):
        return False
    else:
        return os.stat(file_name).st_mtime

# globals should be in the beginning! FIXME
global mdl_latest_time
global sbase_transfer_latest_time
mdl_latest_time = get_file_latest_time(prodrg_xyzin)
sbase_transfer_latest_time = get_file_latest_time(sbase_to_coot_tlc)
# FIXME: this is not a proper name
def mdl_update_timeout_func():

    import operator
    global mdl_latest_time
    global sbase_transfer_latest_time
    
    mdl_now_time   = get_file_latest_time(prodrg_xyzin)
    sbase_now_time = get_file_latest_time(sbase_to_coot_tlc)

    # print "sbase_now_time %s    sbase_latest_time %s" %(sbase_now_time, sbase_transfer_latest_time)

    if (operator.isNumberType(mdl_now_time)):
        if operator.isNumberType(mdl_latest_time):
            if (mdl_now_time > mdl_latest_time):
                mdl_latest_time = mdl_now_time
                import_from_prodrg('mini-prep')

    if operator.isNumberType(sbase_transfer_latest_time):
        if operator.isNumberType(sbase_now_time):
            if (sbase_now_time > sbase_transfer_latest_time):
                sbase_transfer_latest_time = sbase_now_time
                try: # sort of check if file exists?
                    fin = open(sbase_to_coot_tlc, 'r')
                    tlc_symbol = fin.readline()  # need to read more? FIXME
                    fin.close()
                    imol = coot.get_ccp4srs_monomer_and_dictionary(tlc_symbol)
                    if not coot_utils.valid_model_molecule_qm(imol):
                        print "failed to get SBase molecule for", tlc_symbol
                    else:
                        # it was read OK, do an overlap
                        coot_utils.using_active_atom(overlap_ligands, imol,
                                          "aa_imol", "aa_chain_id", "aa_res_no")
                except:
                    print "BL ERROR:: reading sbase file", sbase_to_coot_tlc
                
    return True # return value, keep running; FIXME:: how to stop?


# This is done internally now, by passing lbg a function that runs
# prodrg after the mdl file has been written.  On demand.  We no
# longer need to have a timeout looking for a new prodrg-in.mdl file.
#
# gobject.timeout_add(500, mdl_update_timeout_func)

# return False (if fail) or a list of: the molecule number of the
# selected residue, the prodrg output mol file_name, the prodrg
# output pdb file_name
#
def prodrg_flat(imol_in, chain_id_in, res_no_in):
    """return False (if fail) or a list of: the molecule number of the
    selected residue, the prodrg output mol file_name, the prodrg
    output pdb file_name

    Keyword arguments:
    imol_in -- the molecule number
    chain_id_in -- chain id of the molecule
    res_no_in -- residue number ot the molecule
    
    """

    import operator
    import random

    # return a text string, or at least "" if we can't find the prodrg
    # error message output line.
    def get_prodrg_error_message(log_file):
        if not os.path.isfile(log_file):
            return ""
        else:
            fin = open(log_file, 'r')
            lines = fin.readlines()
            fin.close()
            for line in lines:
                if (" PRODRG: " in line):
                    # assuming only one error?!
                    return line
            return ""  # no error

    selection_string = "//" + chain_id_in + "/" + str(res_no_in)
    imol = coot.new_molecule_by_atom_selection(imol_in, selection_string)
    prodrg_input_file_name = os.path.join("coot-ccp4", "tmp-residue-for-prodrg.pdb")
    prodrg_output_mol_file = os.path.join("coot-ccp4", ".coot-to-lbg-mol")
    prodrg_output_pdb_file = os.path.join("coot-ccp4", ".coot-to-lbg-pdb")
    prodrg_output_lib_file = os.path.join("coot-ccp4", ".coot-to-lbg-lib")
    prodrg_log             = os.path.join("coot-ccp4", "tmp-prodrg-flat.log")
    coot.set_mol_displayed(imol, 0)
    set_mol_active   (imol, 0)
    coot.write_pdb_file(imol, prodrg_input_file_name)
    coot.make_directory_maybe("coot-ccp4")
    arg_list = ["XYZIN",  prodrg_input_file_name,
                "MOLOUT", prodrg_output_mol_file,
                "XYZOUT", prodrg_output_pdb_file,
                "LIBOUT", prodrg_output_lib_file]
    print "arg_list", arg_list
    status = coot_utils.popen_command(cprodrg,
                           arg_list,
                           ["COORDS BOTH", "MINI FLAT", "END"],
                           prodrg_log, True)
    # Does this make sense in python? status being number that is.
    # Maybe rather check for cprodrg exe. FIXME
    if not operator.isNumberType(status):
        coot.info_dialog("Ooops: cprodrg not found?")
        return False
    else:
        if not (status == 0):
            # only for python >=2.5
            mess = "Something went wrong running cprodrg\n" + \
                   get_prodrg_error_message(prodrg_log)
            coot.info_dialog(mess)
            return False
        else:
            # normal return value (hopefully)
            return [imol,
                    prodrg_output_mol_file,
                    prodrg_output_pdb_file,
                    prodrg_output_lib_file]


def prodrg_plain(mode, imol_in, chain_id_in, res_no_in):

    selection_string = "//" + chain_id_in + "/" + \
                       str(res_no_in)
    imol = coot.new_molecule_by_atom_selection(imol_in, selection_string)
    stub = os.path.join("coot-ccp4", "prodrg-tmp-" + str(os.getpid()))
    prodrg_xyzin  = stub + "-xyzin.pdb"
    prodrg_xyzout = stub + "-xyzout.pdb"
    prodrg_cif    = stub + "-dict.cif"
    prodrg_log    = stub + ".log"

    coot.write_pdb_file(imol, prodrg_xyzin)
    result = coot_utils.popen_command(cprodrg,
                           ["XYZIN",  prodrg_xyzin,
                            "XYZOUT", prodrg_xyzout,
                            "LIBOUT", prodrg_cif],
                           ["MINI PREP", "END"],
                           prodrg_log, True)
    coot.close_molecule(imol)
    return [result, prodrg_xyzout, prodrg_cif]


# Why do we pass imol etc and then use active atom?
def fle_view(imol, chain_id, res_no, ins_code):

    import operator
    global have_mingw
    
    # not using active atom, but make property list
    r_flat  = prodrg_flat (imol, chain_id, res_no)
    r_plain = prodrg_plain('mini-no', imol, chain_id, res_no)
    if (r_flat and (r_plain[0] == 0 )):
        imol_ligand_fragment = r_flat[0]
        prodrg_output_flat_mol_file_name = r_flat[1]
        prodrg_output_flat_pdb_file_name = r_flat[2]
        prodrg_output_cif_file_name      = r_flat[3]
        prodrg_output_3d_pdb_file_name   = r_plain[1]
        # 'using_active_atom'
        active_atom = active_residue()
        aa_imol     = active_atom[0]
        aa_chain_id = active_atom[1]
        aa_res_no   = active_atom[2]
        coot.fle_view_internal(aa_imol, aa_chain_id, aa_res_no, "",  # should be from active_atom!!     coot_utils.using_active_atom([[]])
                          imol_ligand_fragment,
                          prodrg_output_flat_mol_file_name,
                          prodrg_output_flat_pdb_file_name,
                          prodrg_output_3d_pdb_file_name,
                          prodrg_output_cif_file_name)
        # touch on Windows!?
        # either distribute touch.exe or use DOS:
        # for non existing file use:
        # copy NUL YourFile.txt
        # for existing:
        # copy /b filename.ext +,,
        #
        # all not neede any more...
        # BL says:: just keep in case I need to touch anything again!!!
        #
        #if (coot_utils.is_windows()):
        #    import subprocess
        #    lbg_ready = os.path.abspath(os.path.join("coot-ccp4",
        #                                              ".coot-to-lbg-mol-ready"))
        #    print "BL DEBUG:: lbg_ready is", lbg_ready
        #    if os.path.isfile(lbg_ready):
        #        subprocess.call("copy /b /y " + lbg_ready + " +,, " + lbg_ready
        #                        , shell=True)
        #    else:
        #        subprocess.call("copy NUL " + lbg_ready, shell=True)
        #else:
        #    coot_utils.popen_command("touch",
        #                  [os.path.join("coot-ccp4",
        #                                ".coot-to-lbg-mol-ready")],
        #                  [],
        #                  "/dev/null", False)

# using cprodrg
# again, why passing imol etc? FIXME
#
def fle_view_to_png(imol, chain_id, res_no, ins_code, neighb_radius,
                    png_file_name):

    # not using active atom, but make property list
    # 'using_active_atom'
    active_atom = active_residue()
    imol     = active_atom[0]
    chain_id = active_atom[1]
    res_no   = active_atom[2]
    
    r_flat  = prodrg_flat (imol, chain_id, res_no)
    r_plain = prodrg_plain('mini-no', imol, chain_id, res_no)
    
    if (r_flat and (r_plain[0] == 0 )):
        imol_ligand_fragment = r_flat[0]
        prodrg_output_flat_mol_file_name = r_flat[1]
        prodrg_output_flat_pdb_file_name = r_flat[2]
        prodrg_output_cif_file_name      = r_flat[3]
        prodrg_output_3d_pdb_file_name   = r_plain[1]
        coot.fle_view_internal_to_png(imol, chain_id, res_no, "",  # should be from active_atom!!     coot_utils.using_active_atom([[]])
                                 imol_ligand_fragment,
                                 prodrg_output_flat_mol_file_name,
                                 prodrg_output_flat_pdb_file_name,
                                 prodrg_output_3d_pdb_file_name,
                                 prodrg_output_cif_file_name, 1,
                                 png_file_name)
                

# import from SRS, callback using sbase_import_function
#
def get_ccp4srs_monomer_and_overlay(comp_id):

    """
    import from SBASE, callback using sbase_import_function
    """

    if (active_residue()):
        with coot_utils.UsingActiveAtom() as [aa_imol, aa_chain_id, aa_res_no,
                                   aa_ins_code, aa_atom_name, aa_alt_conf]:
            imol = coot.get_ccp4srs_monomer_and_dictionary(comp_id)
            overlap_ligands(imol, aa_imol, aa_chain_id, aa_res_no)
            # BL:: again? I dont think so
        coot.get_ccp4srs_monomer_and_dictionary(comp_id)