Warren,
Thanks for the information. So, instead of utilizing the get_dihedral command
I've gone ahead and changed my mode of attack. By using the get_model command
(which is extremely fast), I can access the coordinate information directly.
Then, I wrote a new subroutine that will read in these coordinates and
calculates the necessary dihedral angles from coordinates and NOT from
selections. This drastically improved the performance and I was able to
analyze the sugar pucker in a blink of an eye (as expected). For those who are
interested, the modified code (though lengthier) is below.
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={}
theta={}
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) == 15):
#Construct planes
get_dihedrals(residue,theta)
#Calculate pucker
info = pseudo(residue,theta)
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={}
dihedrals={}
theta={}
#Store new value
store_atom(atom,residue)
else:
store_atom(atom,residue)
first=0
old=new
oldchain=newchain
#Final Residue
#Calculate dihedrals for final residue
if (len(residue) == 15):
#Construct planes
get_dihedrals(residue,theta)
#Calculate pucker for final residue
info = pseudo(residue,theta)
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,theta):
other=2*(sin(math.radians(36.0))+sin(math.radians(72.0)))
#phase=atan2((theta4+theta1)-(theta3+theta0),2*theta2*(sin(math.radians(36.0))+sin(math.radians(72.0))))
phase=atan2((theta['4']+theta['1'])-(theta['3']+theta['0']),theta['2']*other)
amplitude=theta['2']/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):
if (atom.name == "O4'" or atom.name == "O4*"):
residue['O4*X'] = atom.coord[0]
residue['O4*Y'] = atom.coord[1]
residue['O4*Z'] = atom.coord[2]
elif (atom.name == "C1'" or atom.name == "C1*"):
residue['C1*X'] = atom.coord[0]
residue['C1*Y'] = atom.coord[1]
residue['C1*Z'] = atom.coord[2]
elif (atom.name == "C2'" or atom.name == "C2*"):
residue['C2*X'] = atom.coord[0]
residue['C2*Y'] = atom.coord[1]
residue['C2*Z'] = atom.coord[2]
elif (atom.name == "C3'" or atom.name == "C3*"):
residue['C3*X'] = atom.coord[0]
residue['C3*Y'] = atom.coord[1]
residue['C3*Z'] = atom.coord[2]
elif (atom.name == "C4'" or atom.name == "C4*"):
residue['C4*X'] = atom.coord[0]
residue['C4*Y'] = atom.coord[1]
residue['C4*Z'] = atom.coord[2]
return
def get_dihedrals(residue,theta):
C = []
ribose = residue.keys()
ribose.sort()
shift_up(ribose,6)
for i in range (0,12):
C.append(residue[ribose[i]])
theta['0']=dihedral(C)
C = []
shift_down(ribose,3)
for i in range (0,12):
C.append(residue[ribose[i]])
theta['1']=dihedral(C)
C = []
shift_down(ribose,3)
for i in range (0,12):
C.append(residue[ribose[i]])
theta['2']=dihedral(C)
C = []
shift_down(ribose,3)
for i in range (0,12):
C.append(residue[ribose[i]])
theta['3']=dihedral(C)
C = []
shift_down(ribose,3)
for i in range (0,12):
C.append(residue[ribose[i]])
theta['4']=dihedral(C)
return
def shift_up(list,value):
for i in range (0,value):
list.insert(0,list.pop())
return
def shift_down(list,value):
for i in range (0,value):
list.insert(len(list),list.pop(0))
return
def dihedral(C):
DX12=C[0]-C[3]
DY12=C[1]-C[4]
DZ12=C[2]-C[5]
DX23=C[3]-C[6]
DY23=C[4]-C[7]
DZ23=C[5]-C[8]
DX43=C[9]-C[6];
DY43=C[10]-C[7];
DZ43=C[11]-C[8];
#Cross product to get normal
PX1=DY12*DZ23-DY23*DZ12;
PY1=DZ12*DX23-DZ23*DX12;
PZ1=DX12*DY23-DX23*DY12;
NP1=sqrt(PX1*PX1+PY1*PY1+PZ1*PZ1);
PX1=PX1/NP1
PY1=PY1/NP1
PZ1=PZ1/NP1
PX2=DY43*DZ23-DY23*DZ43;
PY2=DZ43*DX23-DZ23*DX43;
PZ2=DX43*DY23-DX23*DY43;
NP2=sqrt(PX2*PX2+PY2*PY2+PZ2*PZ2);
PX2=PX2/NP2
PY2=PY2/NP2
PZ2=PZ2/NP2
DP12=PX1*PX2+PY1*PY2+PZ1*PZ2
TS=1.0-DP12*DP12
if (TS < 0):
TS=0
else:
TS=sqrt(TS)
ANGLE=math.pi/2.0-atan2(DP12,TS)
PX3=PY1*PZ2-PY2*PZ1
PY3=PZ1*PX2-PZ2*PX1
PZ3=PX1*PY2-PX2*PY1
DP233=PX3*DX23+PY3*DY23+PZ3*DZ23
if (DP233 > 0):
ANGLE=-ANGLE
ANGLE=math.degrees(ANGLE)
if (ANGLE > 180):
ANGLE=ANGLE-360
if (ANGLE < -180):
ANGLE=ANGLE+360
return ANGLE
cmd.extend("pucker",pucker)
> Subject: RE: [PyMOL] get_dihedral bottleneck
> Date: Wed, 1 Apr 2009 10:23:33 -0700
> From: war...@delsci.com
> To: magic...@hotmail.com; pymol-users@lists.sourceforge.net
>
> 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.
>
_________________________________________________________________
Create a cool, new character for your Windows Live™ Messenger.
http://go.microsoft.com/?linkid=9656621
