1
2from __future__ import print_function
3
4from pymol.wizard import Wizard
5from pymol import cmd,editor
6from chempy import io
7from copy import deepcopy
8
9import pymol
10import os
11import traceback
12
13src_sele = "_mutate_sel"
14bump_name = "_bump_check"
15
16obj_name = "mutation"
17frag_name = "_tmp_mut"
18mut_sele = "_tmp_mut_sele"
19tmp_obj2 = "_tmp_obj2"
20tmp_sele1 = "_tmp_sele1"
21tmp_sele2 = "_tmp_sele2"
22tmp_hbonds = "_tmp_hbonds"
23
24default_mode = "current"
25default_rep = "lines"
26default_hyd = 'auto'
27default_dep = 'dep'
28default_n_cap = 'none'
29default_c_cap = 'none'
30
31_rot_type_xref = {
32    'GLUH' : 'GLU',
33    'ASPH' : 'ASP',
34    'ARGN' : 'ARG',
35    'LYSN' : 'LYS',
36    'HIP' : 'HIS',
37    'HID' : 'HIS',
38    'HIP' : 'HIS'
39    }
40
41class Mutagenesis(Wizard):
42
43    count = 0
44    cutoff = 3.5
45
46    def __init__(self,_self=cmd):
47        Wizard.__init__(self,_self)
48        cmd=self.cmd
49
50        if self.cmd.get_movie_length() > 0:
51            raise pymol.wizarding.WizardError('Mutagenesis Wizard cannot be used with Movie')
52
53        cmd.unpick()
54
55        self.stored = pymol.Scratch_Storage()
56        self.space = {'stored': self.stored}
57
58        self.bump_scores = []
59        self.dep = default_dep
60
61        self.ind_library = io.pkl.fromFile(os.environ['PYMOL_DATA']+
62                                           "/chempy/sidechains/sc_bb_ind.pkl")
63        self.load_library()
64        self.status = 0 # 0 no selection, 1 mutagenizing
65        self.bump_check = 1
66        self.auto_center = 1
67        self.error = None
68        self.object_name = None
69        self.modes = [
70            'current'
71            ]
72        self.mode = default_mode
73        self.rep = default_rep
74        self.hyd = default_hyd
75        self.n_cap = default_n_cap
76        self.c_cap = default_c_cap
77        residues = list(self.ind_library.keys())
78        # could extent with additional fragments manually as below
79        residues.extend(['GLY','ALA'])
80        residues.extend(['HID','HIE','HIP'])
81        residues.extend(['ARGN','LYSN','ASPH','GLUH'])
82        residues.sort()
83        res_copy = deepcopy(residues)
84        for a in res_copy:
85            residues.append('NT_'+a)
86            residues.append('CT_'+a)
87        self.modes.extend(residues)
88        self.mode_label={}
89        for a in self.modes:
90            self.mode_label[a] = ""+a
91        self.mode_label['current']="No Mutant"
92
93        self.selection_mode = cmd.get_setting_int("mouse_selection_mode")
94        cmd.set("mouse_selection_mode",1)
95
96        smm = []
97        smm.append([ 2, 'Mutant', '' ])
98        smm.append([ 1, 'No change', 'cmd.get_wizard().set_mode("current")' ])
99#        smm.append([ 1, 'N-Term', [] ])
100#        smm.append([ 1, 'C-Term', [] ])
101        smm.append([ 0, '', '' ])
102        for a in self.modes:
103            if a == 'current':
104                pass
105            elif a[0:3]=='NT_':
106                pass
107#                smm[2][2].append([ 1, self.mode_label[a[3:]], 'cmd.get_wizard().set_mode("'+a+'")'])
108            elif a[0:3]=='CT_':
109                pass
110#                smm[3][2].append([ 1, self.mode_label[a[3:]], 'cmd.get_wizard().set_mode("'+a+'")'])
111            else:
112                smm.append([ 1, self.mode_label[a], 'cmd.get_wizard().set_mode("'+a+'")'])
113
114        # group arg, lys, his, glu, asp
115
116        for lst in [ smm ]: # [ smm, smm[2][2], smm[3][2] ]:
117            for a in 'ARG','LYS','HID','GLU','ASP':
118                ix = 0
119                start = 0
120                stop = 0
121                for b in lst:
122                    if start==0:
123                        if b[1][0:]==a:
124                            start = ix
125                            stop = ix + 1
126                    elif b[1][0:3]==a[0:3] or ( b[1][0:2]==a[0:2] and a[0:2]=='HI' ):
127                        stop = ix + 1
128                    ix = ix + 1
129                if start!=0 and stop!=0:
130                    slice = lst[start:stop]
131                    if a != 'HID':
132                        slice2 = [slice[0] ] + [ [0,'',''] ] + slice[1:]
133                        lst[start:stop] = [ [1, self.mode_label[a] + "... " , slice2 ] ]
134                    else:
135                        slice2 = [ slice[3] ] + [ [0,'',''] ] + slice[0:3]
136                        lst[start:stop] = [ [1, self.mode_label['HIS']+ "... ", slice2 ] ]
137
138        self.menu['mode']=smm
139
140
141        self.reps = [
142            'lines',
143            'sticks',
144            'spheres',
145            'dots'
146            ]
147
148        self.rep_name = {
149            'lines' : "Show Lines",
150            'sticks' : "Show Sticks",
151            'spheres' : "Show Spheres",
152            'dots' : "Show Dots",
153            }
154
155        self.dep_name = {
156            'dep' : "Backbone Depen. Rotamers",
157            'ind' : "Backbone Indep. Rotamers"
158            }
159
160        self.hyd_name = {
161            'auto' : "Hydrogens: Current",
162            'keep' : "Hydrogens: Add & Retain",
163#            'polar' : "Polar Hydrogens",
164            'none'  : "Hydrogens: Remove",
165            }
166        self.hyds = [ 'auto', 'keep', 'none' ]
167
168        self.n_cap_name = {
169            'none' : 'Open',
170            'posi' : 'NH3+',
171            'acet' : 'Acetyl',
172            }
173        self.n_caps = [ 'none', 'posi', 'acet' ]
174
175        self.c_cap_name = {
176           'none' : 'Open',
177           'nega' : 'COO-',
178           'amin' : 'Amine',
179           'nmet' : 'N-methyl',
180            }
181        self.c_caps = [ 'none', 'nega', 'amin', 'nmet' ]
182
183        smm = []
184        smm.append([ 2, 'N-Cap', '' ])
185        for a in self.n_caps:
186            smm.append([ 1, self.n_cap_name[a], 'cmd.get_wizard().set_n_cap("'+a+'")'])
187        self.menu['n_cap']=smm
188
189        smm = []
190        smm.append([ 2, 'C-Cap', '' ])
191        for a in self.c_caps:
192            smm.append([ 1, self.c_cap_name[a], 'cmd.get_wizard().set_c_cap("'+a+'")'])
193        self.menu['c_cap']=smm
194
195        smm = []
196        smm.append([ 2, 'Hydrogens', '' ])
197        for a in self.hyds:
198            smm.append([ 1, self.hyd_name[a], 'cmd.get_wizard().set_hyd("'+a+'")'])
199        self.menu['hyd']=smm
200
201        smm = []
202        smm.append([ 2, 'Representation', '' ])
203        for a in self.reps:
204            smm.append([ 1, self.rep_name[a], 'cmd.get_wizard().set_rep("'+a+'")'])
205        self.menu['rep']=smm
206
207        self.deps = [ 'dep', 'ind' ]
208        smm = []
209        smm.append([ 2, 'Rotamers', '' ])
210        for a in self.deps:
211            smm.append([ 1, self.dep_name[a], 'cmd.get_wizard().set_dep("'+a+'")'])
212        self.menu['dep']=smm
213
214        if 'pk1' in cmd.get_names('selections'):
215            cmd.select(src_sele,"(byres pk1)")
216            cmd.unpick()
217            cmd.enable(src_sele)
218            self.status = 1
219            self.error = None
220            self.do_library()
221            cmd.refresh_wizard()
222
223    def load_library(self):
224        if self.dep == 'dep':
225            if not hasattr(self,'dep_library'):
226                self.dep_library = io.pkl.fromFile(os.environ['PYMOL_DATA']+
227                                           "/chempy/sidechains/sc_bb_dep.pkl")
228
229    def set_mode(self,mode):
230        cmd=self.cmd
231        if mode in self.modes:
232            self.mode = mode
233        if self.status==1:
234            self.do_library()
235        cmd.refresh_wizard()
236
237    def set_rep(self,rep):
238        cmd=self.cmd
239        if rep in self.reps:
240            self.rep=rep
241        cmd.hide("("+obj_name+")")
242        cmd.show('lines',obj_name) # always show lines
243        cmd.show(self.rep,obj_name)
244        cmd.refresh_wizard()
245
246    def set_c_cap(self,c_cap):
247        cmd=self.cmd
248        if c_cap in self.c_caps:
249            self.c_cap=c_cap
250        if self.status==1:
251            self.do_library()
252        cmd.refresh_wizard()
253
254    def set_n_cap(self,n_cap):
255        cmd=self.cmd
256        if n_cap in self.n_caps:
257            self.n_cap=n_cap
258        if self.status==1:
259            self.do_library()
260        cmd.refresh_wizard()
261
262    def set_hyd(self,hyd):
263        cmd=self.cmd
264        if hyd in self.hyds:
265            self.hyd=hyd
266        if self.status==1:
267            self.do_library()
268        cmd.refresh_wizard()
269
270    def set_dep(self,value):
271        cmd=self.cmd
272        if value!=self.dep:
273            self.dep = value
274            self.load_library()
275            if src_sele in cmd.get_names("all"):
276                self.do_library()
277            cmd.refresh_wizard()
278
279    def get_panel(self):
280        cmd=self.cmd
281        if int(cmd.get("mouse_selection_mode")!=1):
282            cmd.set("mouse_selection_mode",1)
283        if self.mode == 'current':
284            label = 'No Mutation'
285        else:
286            label = 'Mutate to '+self.mode_label[self.mode]
287        return [
288            [ 1, 'Mutagenesis',''],
289            [ 3, label,'mode'],
290            [ 3, 'N-Cap: '+self.n_cap_name[self.n_cap],'n_cap'],
291            [ 3, 'C-Cap: '+self.c_cap_name[self.c_cap],'c_cap'],
292            [ 3, self.hyd_name[self.hyd],'hyd'],
293            [ 3, self.rep_name[self.rep],'rep'],
294            [ 3, self.dep_name[self.dep],'dep'],
295            [ 2, 'Apply' , 'cmd.get_wizard().apply()'],
296            [ 2, 'Clear' , 'cmd.get_wizard().clear()'],
297            [ 2, 'Done','cmd.set_wizard()'],
298            ]
299
300    def get_event_mask(self):
301        return Wizard.event_mask_pick + Wizard.event_mask_select + Wizard.event_mask_state
302
303    def cleanup(self):
304        cmd=self.cmd
305        global default_mode,default_rep,default_dep,default_hyd
306        global default_n_cap, default_c_cap
307        default_mode = self.mode
308        default_rep = self.rep
309        default_dep = self.dep
310        default_hyd = self.hyd
311        default_n_cap = self.n_cap
312        default_c_cap = self.c_cap
313        cmd.set("mouse_selection_mode",self.selection_mode) # restore selection mode
314        self.clear()
315
316    def clear(self):
317        cmd=self.cmd
318        self.status=0
319        self.bump_scores = []
320        cmd.delete(tmp_hbonds)
321        cmd.delete(tmp_obj2)
322        cmd.delete(mut_sele)
323        cmd.delete(src_sele)
324        cmd.delete(obj_name)
325        cmd.delete(bump_name)
326        cmd.delete("_seeker_hilight")
327        cmd.refresh_wizard()
328
329    def apply(self):
330        cmd=self.cmd
331        if self.status==1:
332            # find the name of the object which contains the selection
333            src_frame = cmd.get_state()
334            try:
335                new_name = cmd.get_object_list(src_sele)[0]
336            except IndexError:
337                print(" Mutagenesis: object not found.")
338                return
339
340            if True:
341                auto_zoom = cmd.get_setting_text('auto_zoom')
342                cmd.set('auto_zoom',"0",quiet=1)
343                if self.lib_mode!="current":
344                    # create copy with mutant in correct frame
345                    state = cmd.get_object_state(new_name)
346                    cmd.create(tmp_obj2, obj_name, src_frame, state)
347                    cmd.set_title(tmp_obj2, state, '')
348                    cmd.color(self.stored.identifiers[4], "?%s & elem C" % tmp_obj2)
349                    cmd.alter(tmp_obj2, 'ID = -1')
350
351                    # select backbone connection atoms
352                    cmd.select(tmp_sele1, 'neighbor ?%s' % (src_sele), 0)
353
354                    # remove residue and neighboring c-cap/n-cap (if any)
355                    cmd.remove("?%s | byres (?%s & "
356                            "(name N & resn NME+NHH | name C & resn ACE))" % (src_sele, tmp_sele1))
357
358                    # create the merged molecule
359                    cmd.create(new_name, "?%s | ?%s" % (new_name, tmp_obj2), state, state)
360
361                    # now connect them
362                    cmd.select(tmp_sele2, '/%s/%s/%s/%s' % ((new_name,) + self.stored.identifiers[:3]))
363                    cmd.bond('?%s & name C' % (tmp_sele1), '?%s & name N' % (tmp_sele2), quiet=1)
364                    cmd.bond('?%s & name N' % (tmp_sele1), '?%s & name C' % (tmp_sele2), quiet=1)
365                    cmd.set_geometry('(?%s | ?%s) & name C+N' % (tmp_sele1, tmp_sele2), 3, 3) # make amide planer
366
367                    # fix N-H hydrogen position (if any exists)
368                    cmd.h_fix('?%s & name N' % (tmp_sele2))
369
370                    # delete temporary objects/selections
371                    cmd.delete(tmp_sele1)
372                    cmd.delete(tmp_sele2)
373                    cmd.delete(tmp_obj2)
374                    self.clear()
375                    # and return to frame 1
376                    cmd.frame(1)
377                    cmd.refresh_wizard()
378                else:
379                    # create copy with conformation in correct state
380                    cmd.create(tmp_obj2,obj_name,src_frame,1)
381
382                    # remove existing c-cap in copy (if any)
383                    cmd.remove("byres (name N and (%s in (neighbor %s)) and resn NME+NHH)"%
384                                (new_name,src_sele))
385                    cmd.remove("(%s) and name OXT"%src_sele)
386
387                    # remove existing n-cap in copy (if any)
388                    cmd.remove("byres (name C and (%s in (neighbor %s)) and resn ACE)"%
389                                (new_name,src_sele))
390
391                    # save existing conformation on undo stack
392#               cmd.edit("((%s in %s) and name ca)"%(new_name,src_sele))
393                    cmd.push_undo("("+src_sele+")")
394                    # modify the conformation
395                    cmd.update(new_name,tmp_obj2)
396#               cmd.unpick()
397                    cmd.delete(tmp_obj2)
398                    self.clear()
399                    # and return to frame 1
400                    cmd.frame(1)
401                    cmd.refresh_wizard()
402                cmd.set('auto_zoom',auto_zoom,quiet=1)
403
404    def get_prompt(self):
405        self.prompt = None
406        if self.status==0:
407            self.prompt = [ 'Pick a residue...']
408        elif self.status==1:
409            self.prompt = [ 'Select a rotamer for %s or pick a new residue...'%self.res_text ]
410        return self.prompt
411
412
413    def do_library(self):
414        cmd=self.cmd
415        pymol=cmd._pymol
416        if not ((cmd.count_atoms("(%s) and name N"%src_sele)==1) and
417                (cmd.count_atoms("(%s) and name C"%src_sele)==1) and
418                (cmd.count_atoms("(%s) and name O"%src_sele)==1)):
419            self.clear()
420            return 1
421        cmd.feedback("push")
422        cmd.feedback("disable","selector","everythin")
423        cmd.feedback("disable","editor","actions")
424        self.prompt = [ 'Loading rotamers...']
425        self.bump_scores = []
426        state_best = 0
427
428        pymol.stored.name = 'residue'
429        cmd.iterate("first (%s)"%src_sele,'stored.name=model+"/"+segi+"/"+chain+"/"+resn+"`"+resi')
430        self.res_text = pymol.stored.name
431        cmd.select("_seeker_hilight",src_sele)
432
433        auto_zoom = cmd.get_setting_text('auto_zoom')
434        cmd.set('auto_zoom',"0",quiet=1)
435        cmd.frame(0)
436        cmd.delete(frag_name)
437        if self.auto_center:
438            cmd.center(src_sele,animate=-1)
439
440        self.lib_mode = self.mode
441        if self.lib_mode=="current":
442            pymol.stored.resn=""
443            cmd.iterate("(%s & name CA)"%src_sele,"stored.resn=resn")
444            rot_type = _rot_type_xref.get(pymol.stored.resn,pymol.stored.resn)
445            if (self.c_cap!='none') or (self.n_cap!='none') or (self.hyd != 'auto'):
446                self.lib_mode = rot_type # force fragment-based load
447            else:
448                cmd.create(frag_name,src_sele,1,1)
449                if self.c_cap=='open':
450                    cmd.remove("%s and name OXT"%frag_name)
451
452        if self.lib_mode!='current':
453            rot_type = self.lib_mode
454            frag_type = self.lib_mode
455            if (self.n_cap == 'posi') and (frag_type[0:3]!='NT_'):
456                if not ( cmd.count_atoms(
457                    "elem C & !(%s) & (bto. (name N & (%s))) &! resn ACE"%
458                                     (src_sele,src_sele))):
459                    # use N-terminal fragment
460                    frag_type ="NT_"+frag_type
461            if (self.c_cap == 'nega') and (frag_type[0:3]!='CT_'):
462                if not ( cmd.count_atoms("elem N & !(%s) & (bto. (name C & (%s))) & !resn NME+NHH"%
463                                     (src_sele,src_sele))):
464                    # use C-terminal fragment
465                    frag_type ="CT_"+frag_type
466            if rot_type[0:3] in [ 'NT_', 'CT_' ]:
467                rot_type = rot_type[3:]
468            rot_type = _rot_type_xref.get(rot_type, rot_type)
469            cmd.fragment(frag_type.lower(), frag_name, origin=0)
470            # trim off hydrogens
471            if (self.hyd == 'none'):
472                cmd.remove("("+frag_name+" and hydro)")
473            elif (self.hyd == 'auto'):
474                if cmd.count_atoms("("+src_sele+") and hydro")==0:
475                    cmd.remove("("+frag_name+" and hydro)")
476            # copy identifying information
477            cmd.alter("?%s & name CA" % src_sele, "stored.identifiers = (segi, chain, resi, ss, color)", space=self.space)
478            cmd.alter("?%s" % frag_name, "(segi, chain, resi, ss) = stored.identifiers[:4]", space=self.space)
479            # move the fragment
480            if ((cmd.count_atoms("(%s & name CB)"%frag_name)==1) and
481                 (cmd.count_atoms("(%s & name CB)"%src_sele)==1)):
482                cmd.pair_fit("(%s & name CA)"%frag_name,
483                             "(%s & name CA)"%src_sele,
484                             "(%s & name CB)"%frag_name,
485                             "(%s & name CB)"%src_sele,
486                             "(%s & name C)"%frag_name,
487                             "(%s & name C)"%src_sele,
488                             "(%s & name N)"%frag_name,
489                             "(%s & name N)"%src_sele)
490            else:
491                cmd.pair_fit("(%s & name CA)"%frag_name,
492                             "(%s & name CA)"%src_sele,
493                             "(%s & name C)"%frag_name,
494                             "(%s & name C)"%src_sele,
495                             "(%s & name N)"%frag_name,
496                             "(%s & name N)"%src_sele)
497
498            # fix the carbonyl position...
499            cmd.iterate_state(1,"(%s & name O)"%src_sele,"stored.list=[x,y,z]")
500            cmd.alter_state(1,"(%s & name O)"%frag_name,"(x,y,z)=stored.list")
501            if cmd.count_atoms("(%s & name OXT)"%src_sele):
502                cmd.iterate_state(1,"(%s & name OXT)"%src_sele,"stored.list=[x,y,z]")
503                cmd.alter_state(1,"(%s & name OXT)"%frag_name,"(x,y,z)=stored.list")
504            elif cmd.count_atoms("(%s & name OXT)"%frag_name): # place OXT if no template exists
505                angle = cmd.get_dihedral("(%s & name N)"%frag_name,
506                                         "(%s & name CA)"%frag_name,
507                                         "(%s & name C)"%frag_name,
508                                         "(%s & name O)"%frag_name)
509                cmd.protect("(%s & name O)"%frag_name)
510                cmd.set_dihedral("(%s & name N)"%frag_name,
511                                 "(%s & name CA)"%frag_name,
512                                 "(%s & name C)"%frag_name,
513                                 "(%s & name OXT)"%frag_name,180.0+angle)
514                cmd.deprotect(frag_name)
515
516
517            # fix the hydrogen position (if any)
518            if cmd.count_atoms("(hydro and bound_to (name N & (%s)))"%frag_name)==1:
519                if cmd.count_atoms("(hydro and bound_to (name N & (%s)))"%src_sele)==1:
520                    cmd.iterate_state(1,"(hydro and bound_to (name N & (%s)))"%src_sele,
521                                      "stored.list=[x,y,z]")
522                    cmd.alter_state(1,"(hydro and bound_to (name N & (%s)))"%frag_name,
523                                    "(x,y,z)=stored.list")
524                elif cmd.select(tmp_sele1,"(name C & bound_to (%s and elem N))"%src_sele)==1:
525                    # position hydro based on location of the carbonyl
526                    angle = cmd.get_dihedral("(%s & name C)"%frag_name,
527                                             "(%s & name CA)"%frag_name,
528                                             "(%s & name N)"%frag_name,
529                                             tmp_sele1)
530                    cmd.set_dihedral("(%s & name C)"%frag_name,
531                                     "(%s & name CA)"%frag_name,
532                                     "(%s & name N)"%frag_name,
533                                     "(%s & name H)"%frag_name,180.0+angle)
534                    cmd.delete(tmp_sele1)
535
536            # add c-cap (if appropriate)
537            if self.c_cap in [ 'amin', 'nmet' ]:
538                if not cmd.count_atoms("elem N & !(%s) & (bto. (name C & (%s))) & !resn NME+NHH"%
539                                       (src_sele,src_sele)):
540                    if cmd.count_atoms("name C & (%s)"%(frag_name))==1:
541                        if self.c_cap == 'amin':
542                            editor.attach_amino_acid("name C & (%s)"%(frag_name), 'nhh', _self=cmd)
543                        elif self.c_cap == 'nmet':
544                            editor.attach_amino_acid("name C & (%s)"%(frag_name), 'nme', _self=cmd)
545                        if cmd.count_atoms("hydro & bound_to (name N & bound_to (name C & (%s)))"%frag_name):
546                            cmd.h_fix("name N & bound_to (name C & (%s))"%frag_name)
547                        # trim hydrogens
548                        if (self.hyd == 'none'):
549                            cmd.remove("("+frag_name+" and hydro)")
550                        elif (self.hyd == 'auto'):
551                            if cmd.count_atoms("("+src_sele+") and hydro")==0:
552                                cmd.remove("("+frag_name+" and hydro)")
553
554            # add n-cap (if appropriate)
555            if self.n_cap in [ 'acet' ]:
556                if not cmd.count_atoms("elem C & !(%s) & (bto. (name N & (%s))) & !resn ACE "%
557                                       (src_sele,src_sele)):
558                    if cmd.count_atoms("name N & (%s)"%(frag_name))==1:
559                        if self.n_cap == 'acet':
560                            editor.attach_amino_acid("name N & (%s)"%(frag_name), 'ace', _self=cmd)
561                        if cmd.count_atoms("hydro & bound_to (name N & bound_to (name C & (%s)))"%frag_name):
562                            cmd.h_fix("name N & (%s)"%frag_name)
563                        # trim hydrogens
564                        if (self.hyd == 'none'):
565                            cmd.remove("("+frag_name+" and hydro)")
566                        elif (self.hyd == 'auto'):
567                            if cmd.count_atoms("("+src_sele+") and hydro")==0:
568                                cmd.remove("("+frag_name+" and hydro)")
569
570
571
572
573        cartoon = (cmd.count_atoms("(%s & name CA & rep cartoon)"%src_sele)>0)
574        sticks = (cmd.count_atoms("(%s & name CA & rep sticks)"%src_sele)>0)
575
576        cmd.delete(obj_name)
577        key = rot_type
578        lib = None
579        if self.dep == 'dep':
580            try:
581                result = cmd.phi_psi("%s"%src_sele)
582                if len(result)==1:
583                    (phi,psi) = list(result.values())[0]
584                    (phi,psi) = (int(10*round(phi/10)),int(10*(round(psi/10))))
585                    key = (rot_type,phi,psi)
586                    if key not in self.dep_library:
587                        (phi,psi) = (int(20*round(phi/20)),int(20*(round(psi/20))))
588                        key = (rot_type,phi,psi)
589                        if key not in self.dep_library:
590                            (phi,psi) = (int(60*round(phi/60)),int(60*(round(psi/60))))
591                            key = (rot_type,phi,psi)
592                    lib = self.dep_library.get(key,None)
593            except:
594                pass
595        if lib is None:
596            key = rot_type
597            lib = self.ind_library.get(key,None)
598            if (lib is not None) and self.dep == 'dep':
599                print(' Mutagenesis: no phi/psi, using backbone-independent rotamers.')
600        if lib is not None:
601            state = 1
602            for a in lib:
603                cmd.create(obj_name,frag_name,1,state)
604                if state == 1:
605                    cmd.select(mut_sele,"(byres (%s like %s))"%(obj_name,src_sele))
606                if rot_type=='PRO':
607                    cmd.unbond("(%s & name N)"%mut_sele,"(%s & name CD)"%mut_sele)
608                for b in a.keys():
609                    if b!='FREQ':
610                        cmd.set_dihedral("(%s & n;%s)"%(mut_sele,b[0]),
611                                         "(%s & n;%s)"%(mut_sele,b[1]),
612                                         "(%s & n;%s)"%(mut_sele,b[2]),
613                                         "(%s & n;%s)"%(mut_sele,b[3]),
614                                         a[b],state=state)
615                    else:
616                        cmd.set_title(obj_name,state,"%1.1f%%"%(a[b]*100))
617                if rot_type=='PRO':
618                    cmd.bond("(%s & name N)"%mut_sele,"(%s & name CD)"%mut_sele)
619                state = state + 1
620            cmd.delete(frag_name)
621            print(" Mutagenesis: %d rotamers loaded."%len(lib))
622            if self.bump_check:
623                cmd.delete(bump_name)
624                cmd.create(bump_name,
625                "(((byobj %s) within 6 of (%s and not name N+C+CA+O+H+HA)) and (not (%s)))|(%s)"%
626                           (src_sele,mut_sele,src_sele,mut_sele),singletons=1)
627                cmd.color("gray50",bump_name+" and elem C")
628                cmd.set("seq_view",0,bump_name,quiet=1)
629                cmd.hide("everything",bump_name)
630                if ((cmd.select(tmp_sele1, "(name N & (%s in (neighbor %s)))"%
631                                (bump_name,src_sele)) == 1) and
632                    (cmd.select(tmp_sele2, "(name C & (%s in %s))"%
633                                (bump_name,mut_sele)) == 1)):
634                    cmd.bond(tmp_sele1,tmp_sele2)
635                if ((cmd.select(tmp_sele1,"(name C & (%s in (neighbor %s)))"%
636                                (bump_name,src_sele)) == 1) and
637                    (cmd.select(tmp_sele2,"(name N & (%s in %s))"%
638                                (bump_name,mut_sele)) == 1)):
639                    cmd.bond(tmp_sele1,tmp_sele2)
640                cmd.delete(tmp_sele1)
641                cmd.delete(tmp_sele2)
642
643                cmd.protect("%s and not (%s in (%s and not name N+C+CA+O+H+HA))"%
644                            (bump_name,bump_name,mut_sele))
645                cmd.sculpt_activate(bump_name)
646                cmd.show("cgo",bump_name)
647                # draw the bumps
648                cmd.set("sculpt_vdw_vis_mode",1,bump_name)
649                state = 1
650                score_best = 1e6
651                for a in lib:
652                    score = cmd.sculpt_iterate(bump_name, state, 1)
653                    self.bump_scores.append(score)
654                    if score < score_best:
655                        state_best = state
656                        score_best = score
657                    state = state + 1
658            cmd.delete(mut_sele)
659        else:
660            cmd.create(obj_name,frag_name,1,1)
661            print(" Mutagenesis: no rotamers found in library.")
662        cmd.set("seq_view",0,obj_name,quiet=1)
663        pymol.util.cbaw(obj_name, _self=cmd)
664        cmd.hide("("+obj_name+")")
665        cmd.show(self.rep,obj_name)
666        cmd.show('lines',obj_name) #neighbor  always show lines
667        if cartoon:
668            cmd.show("cartoon",obj_name)
669        if sticks:
670            cmd.show("sticks",obj_name)
671        cmd.set('auto_zoom',auto_zoom,quiet=1)
672        cmd.delete(frag_name)
673        cmd.frame(state_best)
674
675        # this might be redundant if frame(state_best) changed the state
676        self.do_state(state_best)
677
678        cmd.unpick()
679        cmd.feedback("pop")
680
681    def do_state(self,state):
682        cmd=self.cmd
683        if cmd.get("sculpting")=="on":
684            names = cmd.get_names("all_objects")
685            if (bump_name in names) and (obj_name in names):
686                cmd.update(bump_name,obj_name)
687        if self.bump_scores:
688            print(' Rotamer %d/%d, strain=%.2f' % (state,
689                    cmd.count_states(obj_name), self.bump_scores[state - 1]))
690
691        # update hbonds (polar contacts)
692        cmd.delete(tmp_hbonds)
693        cmd.distance(tmp_hbonds,
694                "?{}".format(obj_name),
695                "(byobj ?{}) and not ?{}".format(src_sele, src_sele),
696                mode=2, state1=state, state2=1)
697
698    def do_select(self,selection):
699        print("Selected!")
700        cmd=self.cmd
701        if (obj_name in cmd.get_names()):
702            if cmd.count_atoms("(%s) and (%s)"%(obj_name,selection)):
703                cmd.deselect()
704                return 1
705        if self.status!=0:
706            cmd.delete(obj_name)
707        cmd.select(src_sele,selection)
708        cmd.unpick()
709        cmd.enable(src_sele)
710        self.status = 1
711        self.error = None
712        self.do_library()
713        cmd.delete(selection)
714        cmd.refresh_wizard()
715        cmd.deselect()
716        return 1
717
718    def do_pick(self,bondFlag):
719        print("Picked!")
720        cmd=self.cmd
721        if bondFlag:
722            self.error = "Error: please select an atom, not a bond."
723            print(self.error)
724        else:
725            if self.status!=0:
726                cmd.delete(obj_name)
727            cmd.select(src_sele,"(byres pk1)")
728            cmd.unpick()
729            cmd.enable(src_sele)
730            self.status = 1
731            self.error = None
732            self.do_library()
733        cmd.refresh_wizard()
734