Re: [PyMOL] get_dihedral bottleneck
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
Re: [PyMOL] get_dihedral bottleneck
On Tue, Mar 31, 2009 at 9:45 PM, Sean Law magic...@hotmail.com wrote: the performance is lightning fast. Can anybody help me improve the overall performance or suggest an alternative? Thanks in advance! A kind of alternative would be to use a custom dihedral command, that avoids evaluating selections: def torsion(a, b, c, d): l = (b - a).cross(c - b) r = (d - c).cross(b - c) a = degrees(l.angle(r)) if (l.cross(r) * (c - b)) 0.0: return -a else: return a where a, b, c, and d are vector classes and 'cross' is defined as usual like: def cross(self, other): return Vector(self.y * other.z - self.z * other.y, self.z * other.x - self.x * other.z, self.x * other.y - self.y * other.x) pymol has a vector class in the chempy package called 'cpv'. gilleain
Re: [PyMOL] get_dihedral bottleneck
= []
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
[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+
