Sean, get_dihedral is expensive because, in the general case, it requires evaluation of four atom selections - a process which is essentially order N with respect to the number of atoms currently loaded inside PyMOL.
As you've discovered, the easiest way to optimize performance is therefore to minimize the number of atoms present in PyMOL by removing all uninvolved atoms. We are working on ways of improving selection performance, at least one of which will be present in 1.2 (the "domain" argument to the select command). Cheers, Warren ________________________________________ From: Sean Law [mailto:magic...@hotmail.com] Sent: Tuesday, March 31, 2009 1:48 PM To: pymol-users@lists.sourceforge.net Subject: [PyMOL] get_dihedral bottleneck Hi PyMol-Community, I have been looking for a PyMOL command to calculate the ribose sugar pucker information (in DNA) but have not found anything. Thus, I began writing a simple Python script that is supposed to take a selection, determine whether or not it contains a ribose sugar ring, and then calculate the phase, amplitude, and pucker orientation for each nucleotide residue within the selection. The completed script works just fine (see the script below) but it slows down significantly in performance when calling cmd.get_dihedral (at least, my crude debugging method using print statements has identified the culprit to be the get_dihedral command) with a LARGE protein-DNA complex (160,000 atoms) loaded. However, if I only load the DNA (with no protein), the performance is lightning fast. Can anybody help me improve the overall performance or suggest an alternative? Thanks in advance! Sean ========================================================= from pymol.cgo import * # get constants from math import * from pymol import cmd def pucker(selection): # Author: Sean Law # Institute: Michigan State University # E-mail: s...@msu.edu obj_name = selection objects = cmd.get_names(selection) for obj in objects: if (obj == selection): obj_name=obj sel=cmd.get_model(selection) first=1 old=" " oldchain=" " residue={} count=0 for atom in sel.atom: new=atom.chain+" "+str(atom.resi) newchain=atom.chain+" "+atom.segi if (not (oldchain == newchain) and first): print " " #Blank line to separate chain output print "%6s %6s %8s Residue" % ("Phase", "Amp", "Pucker") if (not(new == old) and (not first)): #Check that all 5 atoms exist if(len(residue) == 5): #Calculate pucker info = pseudo(residue) print info+" "+old else: print "There is no sugar in this residue" if (not (oldchain == newchain)): print " " #Blank line to separate chain output print "%6s %6s %8s Residue" % ("Phase", "Amp", "Pucker") #Clear values residue={} #Store new value store_atom(atom,residue,obj_name) else: store_atom(atom,residue,obj_name) first=0 old=new oldchain=newchain #Final Residue #Calculate dihedrals for final residue if (len(residue) == 5): #Calculate pucker for final residue info = pseudo(residue) print info+" "+old else: print "There is no sugar in this residue" return def sele_exists(sele): return sele in cmd.get_names("selections"); def pseudo(residue): other=2*(sin(math.radians(36.0))+sin(math.radians(72.0))) theta0=cmd.get_dihedral(residue['C4*'],residue['O4*'],residue['C1*'],residue['C2*']) theta1=cmd.get_dihedral(residue['O4*'],residue['C1*'],residue['C2*'],residue['C3*']) theta2=cmd.get_dihedral(residue['C1*'],residue['C2*'],residue['C3*'],residue['C4*']) theta3=cmd.get_dihedral(residue['C2*'],residue['C3*'],residue['C4*'],residue['O4*']) theta4=cmd.get_dihedral(residue['C3*'],residue['C4*'],residue['O4*'],residue['C1*']) #phase=atan2((theta4+theta1)-(theta3+theta0),2*theta2*(sin(math.radians(36.0))+sin(math.radians(72.0)))) phase=atan2((theta4+theta1)-(theta3+theta0),theta2*other) amplitude=theta2/cos(phase) phase=math.degrees(phase) if (phase < 0): phase=phase+360 # 0 <= Phase < 360 #Determine pucker if (phase < 36): pucker = "C3'-endo" elif (phase < 72): pucker = "C4'-exo" elif (phase <108): pucker = "O4'-endo" elif (phase < 144): pucker = "C1'-exo" elif (phase < 180): pucker = "C2'-endo" elif (phase < 216): pucker = "C3'-exo" elif (phase < 252): pucker = "C4'-endo" elif (phase < 288): pucker = "O4'-exo" elif (phase < 324): pucker = "C1'-endo" elif (phase < 360): pucker = "C2'-exo" else: pucker = "Phase is out of range" info = "%6.2f %6.2f %8s" % (phase, amplitude, pucker) return info def store_atom(atom,residue,obj_name): if (atom.name == "O4'" or atom.name == "O4*"): residue['O4*'] = str(atom_sele(atom,obj_name)) elif (atom.name == "C1'" or atom.name == "C1*"): residue['C1*'] = str(atom_sele(atom,obj_name)) elif (atom.name == "C2'" or atom.name == "C2*"): residue['C2*'] = str(atom_sele(atom,obj_name)) elif (atom.name == "C3'" or atom.name == "C3*"): residue['C3*'] = str(atom_sele(atom,obj_name)) elif (atom.name == "C4'" or atom.name == "C4*"): residue['C4*'] = str(atom_sele(atom,obj_name)) return def atom_sele(atom,obj_name): atom_sele = "" if (not (obj_name == "")): atom_sele = atom_sele+obj_name+" & " if (not (atom.segi == "")): atom_sele = atom_sele+"segi "+atom.segi+" & " if (not (atom.chain == "")): atom_sele = atom_sele+"chain "+atom.chain+" & " if (not (atom.resn == "")): atom_sele = atom_sele+"resn "+atom.resn+" & " if (not (atom.resi == "")): atom_sele = atom_sele+"resi "+str(atom.resi)+" & " if (not (atom.name == "")): atom_sele = atom_sele+"name "+atom.name return atom_sele cmd.extend("pucker",pucker) ________________________________________ Communicate, update and plan on Windows Live Messenger. Get started today.
