Re: [PyMOL] Delete atoms out of unit cell box

[Filipe M. <fi...@xr...>]

I have tweaked a script by Tsjerk Wassenaar that can remove atoms out of =
a=20
unit cell.
I have encountered 2 problems. First is even though i knew the Atom i=20
wanted to remove from the model
i had no easy way to do it because i didn't found any unique identifier=20
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 th=
e=20
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=20
##

def det(a):
     return=20
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=20
##
def minv(a):
     c =3D 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=20
##
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=20
################
##############################################
# p should be the pdb molecule from which=20
the
# unit cell information will be=20
derived

def rmoutside(p =3D None):
     if p =3D=3D None:
         print '# Error: No box specified! #'
         return
     objects =3D 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]) =3D=3D "object:molecule"):
             sel =3D "("+str(objects[j])+")"
             sel =3D cmd.get_model(sel)
             par =3D cmd.get_symmetry(p)
             par[3] =3D par[3]/360*2*math.pi
             par[4] =3D par[4]/360*2*math.pi
             par[5] =3D par[5]/360*2*math.pi

             v_a =3D [par[0],0,0]
             v_b =3D [math.cos(par[5])*par[1],math.sin(par[5])*par[1],0]
             v_c =3D [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 =3D [ v_a,v_b,v_c ]

             # b should now be a box-matrix: [ [ x1, y1, z1 ], [ x2, y2, =
z2=20
], [ x3, y3, z3 ]=20
]
             S =3D [ [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 =3D minv(S)
             list =3D "out_of_box";
             for i in range(len(sel.atom)):
                 sel.atom[i].coord =3D mvmult(T, sel.atom[i].coord)
                 if sel.atom[i].coord[0] < 0 or sel.atom[i].coord[0] > 1 =
or=20
sel.atom[i].coord[1] < 0 or sel.atom[i].coord[1] > 1 or=20
sel.atom[i].coord[2] < 0 or se\
l.atom[i].coord[2] > 1:
                     list =3D "(index "+str(sel.atom[i].index)+" and=20
"+str(objects[j])+") or " + list

                     # this is to try to make sure pymol doesn't seg=20
fault
                     # so when we already have a big list for selection=20
make the selection and empty the list
                     if(len(list) > 2000):
                         cmd.select("out_of_box", list)
                         list =3D "out_of_box"
                     else:
                         sel.atom[i].coord =3D mvmult(S, sel.atom[i].coor=
d)
             cmd.select("out_of_box", list)