I have tweaked a script by Tsjerk Wassenaar that can remove atoms out of a
unit cell.
I have encountered 2 problems. First is even though i knew the Atom i
wanted to remove from the model
i had no easy way to do it because i didn't found any unique identifier
for an Atom in the atom class.
So what i did was to treat each object separately. This seems to work.
The other problem is that pymol segfaults with big selections. I think the
problem is the number of
logical operations in the selection and not the string size.
Here is the script:
from pymol import cmd
import math
## Determining the determinant of matrix a
##
def det(a):
return
a[0][0]*(a[1][1]*a[2][2]-a[2][1]*a[1][2])-a[1][0]*(a[0][1]*a[2][2]-a[2][1]*a[0][2])+a[2][0]*(a[0][1]*a[1][2]-
a[1][1]*a[0][2])
## Basic inversion of a 3x3-matrix a
##
def minv(a):
c = 1/det(a)
return [[ c*(a[1][1]*a[2][2]-a[2][1]*a[1][2]),
-c*(a[0][1]*a[2][2]-a[2][1]*a[0][2]),
c*(a[0][1]*a[1][2]-a[1][1]*a[0][2]) ],
[ -c*(a[1][0]*a[2][2]-a[2][0]*a[1][2]),
c*(a[0][0]*a[2][2]-a[2][0]*a[0][2]),
-c*(a[0][0]*a[1][2]-a[1][0]*a[0][2]) ],
[ c*(a[1][0]*a[2][1]-a[2][0]*a[1][1]),
-c*(a[0][0]*a[2][1]-a[2][0]*a[0][1]),
c*(a[0][0]*a[1][1]-a[1][0]*a[0][1]) ]]
## Multiplying a vector b by matrix a
##
def mvmult(a, b):
return [ a[0][0]*b[0]+a[0][1]*b[1]+a[0][2]*b[2],
a[1][0]*b[0]+a[1][1]*b[1]+a[1][2]*b[2],
a[2][0]*b[0]+a[2][1]*b[1]+a[2][2]*b[2] ]
##############################################
## Protein / image operations
################
##############################################
# p should be the pdb molecule from which
the
# unit cell information will be
derived
def rmoutside(p = None):
if p == None:
print '# Error: No box specified! #'
return
objects = cmd.get_names("objects")
cmd.select("out_of_box","none")
for j in range(len(objects)):
print objects[j]
if(cmd.get_type(objects[j]) == "object:molecule"):
sel = "("+str(objects[j])+")"
sel = cmd.get_model(sel)
par = cmd.get_symmetry(p)
par[3] = par[3]/360*2*math.pi
par[4] = par[4]/360*2*math.pi
par[5] = par[5]/360*2*math.pi
v_a = [par[0],0,0]
v_b = [math.cos(par[5])*par[1],math.sin(par[5])*par[1],0]
v_c = [par[2]*math.cos(par[4]),
par[1]*par[2]*math.cos(par[3])-par[2]*math.cos(par[4])*v_b[0],
math.sqrt(par[2]*par[2]-par[2]*math.cos(par[4])*par[2]*math.cos(par[4])-
par[1]*par[2]*math.cos(par[3])-par[2]*math.cos(par[4])*v_b[0]*par[1]*par[2]*math.cos(par[3])-par[2]*math.cos(par[4])*v_b[0])]
b = [ v_a,v_b,v_c ]
# b should now be a box-matrix: [ [ x1, y1, z1 ], [ x2, y2, z2
], [ x3, y3, z3 ]
]
S = [ [b[0][0], b[1][0], b[2][0]],
[b[0][1], b[1][1], b[2][1]],
[b[0][2], b[1][2], b[2][2]] ]
T = minv(S)
list = "out_of_box";
for i in range(len(sel.atom)):
sel.atom[i].coord = mvmult(T, sel.atom[i].coord)
if sel.atom[i].coord[0] < 0 or sel.atom[i].coord[0] > 1 or
sel.atom[i].coord[1] < 0 or sel.atom[i].coord[1] > 1 or
sel.atom[i].coord[2] < 0 or se\
l.atom[i].coord[2] > 1:
list = "(index "+str(sel.atom[i].index)+" and
"+str(objects[j])+") or " + list
# this is to try to make sure pymol doesn't seg
fault
# so when we already have a big list for selection
make the selection and empty the list
if(len(list) > 2000):
cmd.select("out_of_box", list)
list = "out_of_box"
else:
sel.atom[i].coord = mvmult(S, sel.atom[i].coord)
cmd.select("out_of_box", list)
