On Monday 17 November 2003 14:06, Laurence Pearl wrote:
> I should add that I don't actually know
> how to program in python and have no idea how its 'tuples' work, so
> this is a very FORTRAN-like routine - feel free to improve it.
I have attached a "Pythonized" version of the code. It uses docstrings
for the functions, try...except blocks to check some of the parameters,
and things like ("%8.3f" * 17) to save some repetition of format
specifiers. Python is cool!
I have also changed the script to add two commands to PyMol using
cmd.extend: "camtrav first,nframes,selection,zoom,buffer" and
"seqfly object,firstres,lastres". Usage examples:
#---------------------------------
load protein.pdb
run camera_travel.py
mset 1 x150
zoom all
camtrav 1,100,(resi 50),1,4
mplay
# This will zoom in from the complete molecule to residue 50.
#---------------------------------
#---------------------------------
load protein.pdb
run camera_travel.py
mset 1 x150
seqfly protein,1,8
mplay
# This will visit the sequence from residue 1 to 8
#---------------------------------
Finally, I ran into some trouble in a lysozyme structure (PDB 103L). At
residue 139, cosom becomes slightly larger than 1 (1.0008 or so) and
acos(cosom) fails. Also, omega and sinom become zero in this case, so
the divisions by sinom fail. I added some checks to normalize this
situation and stop the code from crashing, but we should look into
what's causing this.
--
Lieven Buts
Department of Ultrastructure
Vrije Universiteit Brussel# camera_travel - Laurence Pearl, November 2003
# Pythonized by Lieven Buts, November 2003
import cmd
from math import sqrt,acos,sin
def quaternion(view):
"""
Returns a quaternion representation of a view matrix.
"""
nxt = [1,2,0]
q = [0.0,0.0,0.0,1.0]
tr = view[0]+view[4]+view[8]
if tr > 0.0 :
s = sqrt(tr + 1.0)
qw1 = s / 2.0
s = 0.5 / s
return ( (view[5] - view[7]) * s,
(view[6] - view[2]) * s,
(view[1] - view[3]) * s,
qw1 )
else :
i = 0
if (view[4] > view[0]): i = 1
if (view[8] > view[i+3*i]): i = 2
j = nxt[i]
k = nxt[j]
s = sqrt ((view[i+i*3] - (view[j+j*3] + view[k+k*3])) + 1.0)
q[i] = s * 0.5
if (s != 0.0): s = 0.5 / s
q[3] = (view[k+3*j] - view[j+3*k]) * s
q[j] = (view[j+3*i] + view[i+3*j]) * s
q[k] = (view[k+3*i] + view[i+3*k]) * s
return q
def camera_travel(first,nframes=30,sel='(all)',zflag=0,zlevel=2):
"""
Generate progressive view matrices to move the camera smoothly
from the current view to a view defined by a PyMol selection.
first - start frame
nframes - duration
sel - atom selection that defines the orientation at the end of the
sequence
zflag - flag to indicate whether the final view should be 'zoomed'
or not
zlevel - buffer space for zooming final view (as in cmd.zoom)
"""
try:
first=int(first)
except:
print "camera_travel: first frame must be an integer"
return -1
try:
nframes=int(nframes)
except:
print "camera_travel: number of frames"
ff=float(1.0/nframes)
old_view = cmd.get_view(2)
# print ( "view : (" + "%8.3f, "*17 + "%8.3f)" ) % (old_view)
# print "oldtran : %8.3f %8.3f %8.3f" % (old_view[12], old_view[13],
old_view[14])
# do orient operation on selection
cmd.orient(sel,0)
# if zoom to selection is required add this into view matrix
if zflag:
cmd.zoom(sel,zlevel,0,1)
# get new view
new_view = cmd.get_view()
# capture new zoom/clip parameters
ozc1 = new_view[11]
ozc2 = new_view[15]
ozc3 = new_view[16]
# calculate shift in zoom/clip parameters
dzc1 = (ozc1 - old_view[11]) * ff
dzc2 = (ozc2 - old_view[15]) * ff
dzc3 = (ozc3 - old_view[16]) * ff
ozc1 = old_view[11]
ozc2 = old_view[15]
ozc3 = old_view[16]
# capture new translation vector component
ox = new_view[12]
oy = new_view[13]
oz = new_view[14]
# calculate shift vector
dx = ox - old_view[12]
dy = oy - old_view[13]
dz = oz - old_view[14]
dx = dx*ff
dy = dy*ff
dz = dz*ff
ox = old_view[12]
oy = old_view[13]
oz = old_view[14]
# capture old and new rotation matrix components in quaternion form
# m[0][0] = v[0] m[0][1] = v[1] m[0][2] = v[2]
# m[1][0] = v[3] m[1][1] = v[4] m[1][2] = v[5]
# m[2][0] = v[6] m[2][1] = v[7] m[2][2] = v[8]
qx1,qy1,qz1,qw1 = quaternion(old_view)
qx2,qy2,qz2,qw2 = quaternion(new_view)
# calc cosine
cosom = qx1 * qx2 + qy1 * qy2 + qz1 * qz2 + qw1 * qw2
limit = 0.001
if cosom>1.0+limit:
raise ValueError,"Cosine of omega way out of range (positive)"
elif cosom>1.0:
print "Warning: cosom corrected from ",cosom,"to",
cosom = 1.0
print cosom
if cosom<-1.0-limit:
raise ValueError,"Cosine of omega way out of range (negative)"
elif cosom<-1.0:
print "Warning: cosom corrected from ",cosom,"to",
cosom = 1.0
print cosom
# adjust signs
if (cosom < 0.0):
cosom = -cosom
to0 = -qx2
to1 = -qy2
to2 = -qz2
to3 = -qw2
else:
to0 = qx2
to1 = qy2
to2 = qz2
to3 = qw2
# calc coefficients
omega = acos(cosom)
sinom = sin(omega)
if sinom==0.0:
sinom=limit
print "Warning: sinom corrected!"
# restore old view
cmd.set_view( ("%8.3f, " * 17 + "%8.3f") % tuple(old_view) )
# loop interpolating over nframes generating interpolated quaternion
for a in range(nframes+1):
scale0 = sin((1.0 - float(a*ff)) * omega) / sinom
scale1 = sin(float(a*ff) * omega) / sinom
rx = scale0 * qx1 + scale1 * to0;
ry = scale0 * qy1 + scale1 * to1;
rz = scale0 * qz1 + scale1 * to2;
rw = scale0 * qw1 + scale1 * to3;
# convert back to matrix
x2 = rx + rx
y2 = ry + ry
z2 = rz + rz
xx = rx * x2
xy = rx * y2
xz = rx * z2
yy = ry * y2
yz = ry * z2
zz = rz * z2
wx = rw * x2
wy = rw * y2
wz = rw * z2
nv0 = 1.0 - (yy + zz)
nv3 = xy - wz
nv6 = xz + wy
nv1 = xy + wz
nv4 = 1.0 - (xx + zz)
nv7 = yz - wx
nv2 = xz - wy
nv5 = yz + wx
nv8 = 1.0 - (xx + yy)
# update translation vector
ox = ox + dx
oy = oy + dy
oz = oz + dz
# update zoom/clip parameters if required
if zflag:
ozc1 = ozc1 + dzc1
ozc2 = ozc2 + dzc2
ozc3 = ozc3 + dzc3
cmd.mdo("%d" % (first), ("set_view (" + "%8.3f, "*17 +
"%8.3f)") %
(nv0,nv1,nv2,nv3,nv4,nv5,nv6,nv7,nv8,old_view[9],old_view[10],ozc1,
ox,oy,oz,ozc2,ozc3,old_view[17]))
first += 1
cmd.set_view( ("%8.3f, "*17 + "%8.3f") %
(nv0,nv1,nv2,nv3,nv4,nv5,nv6,nv7,nv8,old_view[9],old_view[10],ozc1,
ox,oy,oz,ozc2,ozc3,old_view[17]))
def fly_sequence(object,first,last):
"""
Visit all residues in an object consecutively.
Requires 15 movie frames per residue.
object -
first -
last -
"""
fr=1
for ires in range(int(first),int(last)+1):
res = "(" + object + (" and resi %d" % (ires,)) + ")"
camera_travel(fr,9,res,1,2)
fr = fr + 15
# Install new PyMol commands
cmd.extend("seqfly",fly_sequence)
cmd.extend("camtrav",camera_travel)
