Re: [PyMOL] you know what would be cool?

2003-11-18 Thread Lieven Buts
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] 

Re: [PyMOL] you know what would be cool?

2003-11-18 Thread Laurence Pearl
I'm still not sure about all this 'cool' stuff - Python looks to me 
like the illegitimate spawn of C and BASIC, but then I used to program 
in 6502 machine code so what do I know .


Below is a variant of camera_travel (as Pythonised by Lieven Buts) 
which takes an explicit final view matrix as an argument rather than a 
PyMol selection as in the original version. This allows you to specify 
the exact view you want, rather than relying on the combination of 
PyMol 'orient' and 'zoom' command to generate it for you. This is 
particularly useful for defining a base view to which you can return.


For those with tricky E-mail clients the routine is also attached as a 
text file.


Cheers,

Laurence



def camera_view_travel(first,nframes=30,new_view=(\
 1.0,0.0,0.0,\
 0.0,1.0,0.0,\
 0.0,0.0,1.0,\
 0.0,0.0,  -50.0,\
 0.0,0.0,0.0,\
40.0,  100.0,0.0 )):

Generate progressive view matrices to move the camera smoothly
from the current view to a new view provided as an argument.

first   - start frame
nframes - duration
new_view  - PyMol view matrix that defines the view at the end 
of the sequence



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])


#   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 cosom1.0+limit:
raise ValueError,Cosine of omega way out of range (positive)
elif cosom1.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 *