generate1.py

# pymol -c generate1.py

# NOTE: obsolete -- PyMOL now uses Dunbrack rotamers by default

from chempy import io
from glob import glob

sys.path.append(".")

from pymol import cmd

import pymol
from pymol import cmd

def find_torsions(res_sele):
   # this routine gets a list of atom ids which correspond to the
   # free torsions in the sidechain
   res_sele = "("+res_sele+")"
   if cmd.count_atoms("(%s and n;N,CA)"%res_sele,quiet=1)<2:
      return []
   result = []
   done = 0
   # set up the starting trunk
   cmd.select("trunk","(%s and n;N,CA,C,O)"%res_sele)
   first_id = cmd.identify("(%s and n;N)"%res_sele)[0]
   second_id = cmd.identify("(%s and n;CA)"%res_sele)[0]
   to_do_list = [(first_id,second_id)]
   while len(to_do_list):
      (first_id,second_id) = to_do_list.pop(0)
      # find candidates for the 3rd atom in the dihedral
      cmd.select("candi","(%s and (neighbor trunk))"%res_sele)
      pymol.stored.list = []
      cmd.iterate("candi","stored.list.append((name,ID))")
      candi_list = pymol.stored.list
      if len(candi_list):
         # we may have a dihedral
         dihe_list = []
         for a in candi_list:
            if cmd.select("termi","(%s and (not trunk) and (neighbor id %d))"%
                          (res_sele,a[1]))>0:
               # apparently we do, but is it acyclic
               # (i.e., can edit split the molecule?)
               cmd.edit("(%s and trunk and (neighbor id %d))"%(res_sele,a[1]),
                        "(%s and id %d)"%(res_sele,a[1]))
               if 'pkfrag2' in cmd.get_names('selections'):
                  # yes, so this is a valid dihedral
                  dihe_list.append(a)
               else:
                  # no, so add the third atom to the trunk
                  cmd.select("trunk","(trunk or (%s and id %d))"%(res_sele,a[1]))                  
            else:
               # no fourth atom, so add this atom to the trunk
               cmd.select("trunk","(trunk or (%s and id %d))"%(res_sele,a[1]))
         if len(dihe_list):
            # choose 3rd atom using alphanumeric order
            dihe_list.sort()
            third_id = dihe_list[0][1]
            # if there is another third atom, then repeat later      
            if len(dihe_list)>1: 
               to_do_list.insert(0,(first_id,second_id))
            # now choose the 4th atom, which we know exists, using a similar criterion
            cmd.select("termi","(%s and (not trunk) and (neighbor id %d))"%
                       (res_sele,third_id))
            pymol.stored.list=[]
            cmd.iterate("termi","stored.list.append((name,ID))")
            termi_list = pymol.stored.list
            termi_list.sort()
            fourth_id = termi_list[0][1]
            # at this point, we should have a complete dihedral
            # add the third atom into the trunk, and store the second and third for
            # outward extension later on
            to_do_list.append((second_id,third_id))
            cmd.select("trunk","(trunk or (%s and id %d))"%(res_sele,third_id))      
#            cmd.show('sticks','(%s and (id %d,%d,%d,%d))'%
#                     (res_sele,first_id,second_id,third_id,fourth_id))
            result.append((first_id,second_id,third_id,fourth_id))
   return result

def measure_all_torsions(sele,prot_dict):
   id_list = cmd.identify("((%s) and name ca)"%sele)
   for a in id_list:
      res_sele = "(alt '' and (byres id %d) and (%s))"%(a,sele)
      lst = find_torsions(res_sele)
      if len(lst):
         cmd.iterate("(%s and n;ca)"%res_sele,"stored.resn=resn")
         resn = pymol.stored.resn
         print resn
         name_list = []
         for a in lst:
            pymol.stored.names = []
            for b in a:
               cmd.iterate("(%s and id %d)"%(res_sele,b),"stored.names.append(name)")
            name_list.append(tuple(pymol.stored.names))
         res_dict = {}
         for a in name_list:
            dihe = cmd.get_dihedral("(%s and name %s)"%(res_sele,a[0]),
                                    "(%s and name %s)"%(res_sele,a[1]),
                                    "(%s and name %s)"%(res_sele,a[2]),
                                    "(%s and name %s)"%(res_sele,a[3]))
            res_dict[a]=dihe
         if not prot_dict.has_key(resn):
            prot_dict[resn]=[]
         prot_dict[resn].append(res_dict)

sc_dict = {}

cmd.feedback('disable','selector','result')
cmd.feedback('disable','executive','action')

for fyle in glob("source/*"):
   cmd.load(fyle)
   cmd.remove("(het)")
   cmd.remove("(hydro or not alt '')")
   cmd.remove("(byres b>25)") # only use well-defined residues
   cmd.remove("(not (byres name ca))")
   cmd.remove("(not (byres name n))")
   cmd.remove("(not (byres name c))")   
   cmd.refresh()
   measure_all_torsions("all",sc_dict)
   cmd.delete("all")
   
#for a in sc_dict.keys():
#   print a,dict[a]
#print len(dict.keys())
      
io.pkl.toFile(sc_dict,"sc_raw.pkl")