Pokazywanie postów oznaczonych etykietą Pymol. Pokaż wszystkie posty
Pokazywanie postów oznaczonych etykietą Pymol. Pokaż wszystkie posty

niedziela, maja 15, 2011

My first cover

Two month ago I was asked to design a cover image for a small polish chemical journal - LAB. If I may quote the editor:
Please prepare something attractive, i.e. the molecule should be very colourful and elaborate.
As I like to play with graphics, especially when its chemistry related, I'd accepted the offer. I thought for some time what molecule should I use for the cover, when a colleague whose work is opioid related proposed naloxone. Its a drug used in emergency treatments for opioid overdose (mostly heroin). I've liked it from the first sight so I searched the protein data base for a suitable model to do the visualization. I have found 1MX9 which is a complex of human liver carboxylesterase with a naloxone derivative. Here is my final work that I've sent to the editor: After a month when I went to work a package was lying at my desk. Containing this magazine. The scan below does not reflect the actual appearance because of the print on a semi gloss paper. It looks more darker then it is. Neither the less it feels really good to see my own work in print.

poniedziałek, stycznia 31, 2011

Grouping conformations in pymol

Lately I had to do some visualizations of molecules after MD written in multirecord PDB files, each file containing MD frames for a different molecule. Loading it and preparing the conformations was no big deal - I just split every conformation into different object.
split_states D1_wszystko The command above will generate new objects called D1_wszystko_0001, D1_wszystko_0002, D1_wszystko_0003 and so on. Since the goal was to compare the different structures and because there was over six hundred of them I had to use some python magic to make my life easier. Lucky for me Pymol comes with a grouping command.
for i in range(1,624): cmd.group("D1_AMBER","D1_wszystko_%04d" % i,"add") This will make a group called D1_AMBER containing all the D1_wszystko_XXXX objects that could be easily collapsed when not needed.

środa, grudnia 24, 2008

Human Adenovirus 2

During the holidays I wanted to try something out that was longer on my mind - remote rendering. At the university I have a pretty plain-desktop machine with Archlinux installed on it. It's not much. Integrated Intel graphic card and 64 bit 2GHz CPU. Nothing unusual, but I was curious how do the limits of the machine look if it comes to render some bigger models. Let's say virus capsids. My previous experiments with my home pc on this field failed on more complicated then ball representations. To say it simply - I run out of memory. It doesn't really matters if it's Linux or Windows, because on a 32 bit machine you have only maximum 4 GB of memory available. On a 64 bit system, however there's theoretically 16 EB of memory available so rendering viruses shouldn't be a problem. I went to VIPERdb and picked the first interesting looking virus I found. It was the Human Adenovirus 2. It's a medium-sized (90–100 nm) virus which is responsible for upper respiratory infections in children. I've downloaded both the half and full size capsid models, uploaded them to the remote machine and extracted the content. After that I went first for some low resource rendering gradually increasing the load of computation needed to be done. Here is how a ball model of the HAd2 capsid with b-factor colouring looks like:
It took about 5 minutes to render. By changing the representation mode from spheres to surface the render time increased to about half an hour. Surface representation also requires significant more memory so I added to the systems 3GB (1GB RAM + 2GB of swapspace) another 8GB by creating a swap file.
# cd /
# dd if=/dev/zero of=swapfile bs=1024 count=8388608
# chmod 600 swapfile
# mkswap swapfile
# swapon swapfile
As it showed up later this still wasn't enough. I've ended up with creating a second swapfile during the rendering and a total of 22 GB memory.
#free -m
             total       used       free     shared    buffers     cached
Mem:           994        843        151          0        118        410
-/+ buffers/cache:        314        680
Swap:        21295          2      21292
This is how a half capsid looks after rendering with surface (also b-factor colouring but this time with a different colour scheme):
And here's the full one with a better surface quality:
It took almost 45h to render and consumed ~19GB of memory. So far it was the longest and most resource consuming rendering I have ever done, but I think there's pretty much room for improvement ];) All the render was done by writing a render script and then run in the background on the remote machine by using the unix nohup command. nohup pymol -c render.py & This is necessary to prevent the application from receiving a hangup signal after logout. The "-c" flag turns pymol to command-mode so it wouldn't use any GUI.
Last but not least here's a small animation of the full capsid.
[Update] I'm starting to hate blogger for what he is doing with the quality of the uploaded movies.

wtorek, grudnia 23, 2008

Building Pymol SVN on Archlinux

Lately I've been neglecting my blog for the sake of my work as a Ph.D. student (programming in Fortran). I'm going to catch up on this. Writing some Fortran code isn't that hard even if you are starting to learn the language from scratch like I do, but definitely it isn't the type of "coding is fun" programming language like Python for example. Speaking of Python, since 3rd of December your python scripts are getting old, because Python 3.0 is out. The new version brings major changes in the language that make you previous HelloWorld.py program stop working (print is now a function and should be called out print("Hello World!") instead of print "Hello World!"). Many programs I use require Python 2.4 or 2.5, some like OpenOffice.org Python 2.6 and soon some will require Python 3.0. Keeping compatibility on a system is getting problematic. Luckily I've installed Archlinux on my desktop at work which makes life easier if you know how to use it. I wanted to install the latest Pymol version from svn, but I've found it's incompatible with Python 2.6 which was required as I mentioned above for OO.o. So here's what I've done.

First I've installed a separate python version from the repositories - Python 2.4 # pacman -Sy python24
After this I needed Pmw for the GUI. It can be done by building from source like this:
# wget http://downloads.sourceforge.net/pmw/Pmw.1.3.2.tar.gz
# tar -xvzf Pmw.1.3.2.tar.gz
# cd Pmw.1.3.2/src
# python2.4 setup.py build
# python2.4 setup.py install
of by using a modified python-pmw PKGBUILD from the AUR. Just make sure it uses python 2.4 by changing in the build function the default python:
build() {
cd $startdir/src/Pmw.$pkgver/src
python setup.py build
python setup.py install --prefix=/usr --root=${startdir}/pkg
}
to python2.4:

build() {
cd $startdir/src/Pmw.$pkgver/src
python2.4 setup.py build
python2.4 setup.py install --prefix=/usr --root=${startdir}/pkg
}

After this you are ready to build pymol from the SVN:
# svn co https://pymol.svn.sourceforge.net/svnroot/pymol/trunk/pymol pymol
# cd pymol
# python2.4 setup.py install
# python2.4 setup2.py install
# cp ./pymol /usr/local/bin/pymol

Thats It!
Now an upgrade of Python to 3.0 will not mess up with your Pymol installation.

środa, września 17, 2008

Cute molecules

For me is Pymol the number one tool for molecule visualization. Sure, there's a lot of visualization software out there on the net and some of them are quite impressing, but I haven't found one that is so flexible and extensible like Pymol. Basically almost all chemistry visualizations tools are OpenGL applications that squeeze more or less the last pixel from your GFX-card. It's just a matter of the settings. So hat makes the difference? The answer is Python. The core program is written in C, but above this is a bunch of Python modules and that means the program is hackable without recompiling. All you need is a plain text editor...

Ever heard about QuteMol? It's a very nice program and the images are just astonishing. The reason why I didn't convinced myself to try it is because it's only available for the Windows and Mac platforms and I'm starting to use my Arch Linux more and more. Ok, now you are probably asking yourself what I asked myself - how do I get such images in Pymol? I was not the first who thought about it. Here's a pymolwiki recipe how it can be done. Like I sed earlier - it's just a matter of the settings. You could write a script for it but I've decided to go one step further integrating it into Pymol. After some cups of coffee I've started playing around with the Pymol GUI python modules and found the right ones - menu.py and preset.py. All modifications were done on a linux Pymol version 0.99rc6, but besically it should also work on other Pymol versions.
First I've opened the file preset.py and defined my own representation preset function called qutemollike (add the code below at the end of the file):

preset.py(Toggle Plain Text)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
def qutemollike(selection="(all)",_self=cmd):
    s = tmp_sele
    cmd.set_color("oxygen",[1.0,0.4,0.4])
    cmd.set_color("nitrogen", [0.5,0.5,1.0])
    cmd.remove("solvent")
    cmd.as("spheres")
    util.cbaw()
    cmd.bg_color("white")
    cmd.set("light_count",10)
    cmd.set("spec_count",1)
    cmd.set("shininess",10)
    cmd.set("specular",0.25)
    cmd.set("ambient",0)
    cmd.set("direct",0)
    cmd.set("reflect",1.5)
    cmd.set("ray_shadow_decay_factor",0.1)
    cmd.set("ray_shadow_decay_range", 2)
    cmd.unset("depth_cue")

and then added a callback (line #495) to it in the file menu.py witch describes all the menus used in the viewer window.

menu.py(Toggle Plain Text)

482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
def presets(s):
    return [[ 2, 'Preset:'       ,''                        ],     
              [ 1, 'simple'   ,'preset.simple("'+s+'")'          ],
              [ 1, 'simple (no solvent)'   ,'preset.simple_no_solv("'+s+'")'          ],           
              [ 1, 'ball and stick' , 'preset.ball_and_stick("'+s+'")' ],
              [ 1, 'b factor putty' , 'preset.b_factor_putty("'+s+'")' ],
              [ 1, 'technical'   , 'preset.technical("'+s+'")'          ],
              [ 1, 'ligands'   , 'preset.ligands("'+s+'")'          ],
              [ 1, 'ligand sites'   , preset_ligand_sites(s)         ],
              [ 1, 'pretty ', 'preset.pretty("'+s+'")'          ],
              [ 1, 'pretty (with solvent)'     , 'preset.pretty_solv("'+s+'")'          ],
              [ 1, 'publication '   , 'preset.publication("'+s+'")'          ],
              [ 1, 'publication (with solvent)'   , 'preset.pub_solv("'+s+'")'          ],
              [ 1, 'QuteMol like', 'preset.qutemollike("'+s+'")' ],
              [ 0, ''               ,''                             ],                      
              [ 1, 'default'   ,'preset.default("'+s+'")'          ],           
              ]


Here is a screenshot showing the added "QuteMol like" option in the preset submenu.

As you can see I've loaded a caffeine molecule. I give it the final polish by typing in the console: set ray_trace_mode,1
ray 800,600
And here's how it looks now.

niedziela, czerwca 29, 2008

More "binding site apples"

It took me more time then I've expected, but finally I've updated my previous script. And here it is.

bp_surface.py(Toggle Plain Text)

# -*- coding: utf-8 -*-

# -------------------------------------------------------------------------------
# bp_surface.py - Binding pocket surface generator for Pymol version 2.0
# -------------------------------------------------------------------------------

# Copyright (C) 2008 by Lightnir - lightnir@gmail.com

# This script is free software; you can redistribute it and#or modify
# it under the terms of the GNU Library General Public License as
# published by the Free Software Foundation; either version 2 of the
# License, or (at your option) any later version.

# This script is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.

# You should have received a copy of the GNU Library General Public
# License along with this program; if not, write to the
# Free Software Foundation, Inc.,
# 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

from pymol import cmd
from Tkinter import *
from pymol.cgo import *
from math import *
import Pmw
import tkColorChooser

BPSFrame = None

def __init__(self):
    self.menuBar.addcascademenu('Plugin', 'LightnirsPlugins', 'Lightnirs Plugins',
                                    label='Lightnirs Plugins')
    self.menuBar.addmenuitem('LightnirsPlugins', 'command',
                                 'BPsurface',
                                 label='BP surface',
                                 command = lambda s=self: BPSurface(s))
class BPSurface:
    # Class variables
    selections = []
    cv = IntVar()
    pp = IntVar()
    cp = IntVar()
    rv = IntVar()
    colorize = IntVar()
    Default_Color = None
    
    res_states = [['Hydrophobic','Gray',0], ['Aromatic','White',0], ['Polar','Green',0], ['Positive','Red',0], ['Negative','Blue',0], ['Cysteine','Yellow',0], ['Proline','Magenta',0]]

    def __init__(self,app):
        global BPSFrame

        if BPSFrame==None:
            BPSFrame = Toplevel(width=400, height=300)              # create a root window
            BPSFrame.title("BP surface editor")                     # Set window title
            BPSFrame.resizable(0, 0)                                # Disable window resize
            BPSFrame.geometry('-40+40')                             # Set window placement
            self.get_selections()
            self.Default_Color = BPSFrame.cget("bg")
            # Create Ligand selection Combobox
            BPSFrame.cb1 = Pmw.ComboBox(BPSFrame,
                label_text = 'Ligand:',
                labelpos = 'nw',
                scrolledlist_items = BPSurface.selections
            )
            BPSFrame.cb1.grid(row=0, column=0,columnspan=2)

            # Create binding pocket Combobox
            BPSFrame.cb2 = Pmw.ComboBox(BPSFrame,
                label_text = 'Binding pocket:',
                labelpos = 'nw',
                scrolledlist_items = BPSurface.selections
            )
            BPSFrame.cb2.grid(row=1, column=0,columnspan=2)

            # Create Z axis Slider
            BPSFrame.slider1 = Scale(BPSFrame, from_=-20, to=20, command=self.scaling, tickinterval=10, resolution=0.5, label='Z axis cutoff', orient=HORIZONTAL)
            BPSFrame.slider1.grid(row=2, column=0, columnspan=2, sticky=N+W+S+E,)

            # Create X axis Slider
            BPSFrame.slider2 = Scale(BPSFrame, from_=-90, to=90, command=self.scaling, tickinterval=30, resolution=1, label='X axis rotation',orient=HORIZONTAL)

            # Create Y axis Slider
            BPSFrame.slider3 = Scale(BPSFrame, from_=-90, to=90, command=self.scaling, tickinterval=30, resolution=1, label='Y axis rotation', orient=HORIZONTAL)

            # Add Checkboxes
            BPSurface.cv.set(1)
            cbx = Checkbutton(BPSFrame, text="Create new object", command=None, variable=BPSurface.cv, onvalue="1", offvalue="0")
            cbx.grid(row=5,column=0, columnspan=2, sticky=W)
            self.CreatePlane()
            
            BPSurface.pp.set(1)
            cbx2 = Checkbutton(BPSFrame, text="Plane preview", command=self.TogglePlane, variable=BPSurface.pp, onvalue="1", offvalue="0")
            cbx2.grid(row=6,column=0, columnspan=2, sticky=W)

            BPSurface.cp.set(0)
            cbx3 = Checkbutton(BPSFrame, text="Use custom plane", command=self.ToggleCustomPlane, variable=BPSurface.cp, onvalue="1", offvalue="0")
            cbx3.grid(row=7,column=0, columnspan=2, sticky=W)

            BPSurface.rv.set(0)
            cbx4 = Checkbutton(BPSFrame, text="Make residue selections", command=self.ToggleResidueSelect, variable=BPSurface.rv, onvalue="1", offvalue="0")
            cbx4.grid(row=8,column=0, columnspan=2, sticky=W)

            # Add some buttons
            btn1=Button(BPSFrame, text='Refresh', padx=0, pady=0, command = self.refresh_list)
            btn1.grid(row=9, column=0, sticky=N+W+S+E, padx=0, pady=0,ipadx=0,ipady=0)
            btn2=Button(BPSFrame, text='Clear', padx=0, pady=0, command = lambda r=1,d=1: self.clear_flags(r,d))
            btn2.grid(row=9, column=1, sticky=N+W+S+E, padx=0, pady=0,ipadx=0,ipady=0)
            btn3=Button(BPSFrame, text='Apply', padx=0, pady=0, command = self.ok)
            btn3.grid(row=10, column=0, sticky=N+W+S+E,)
            btn4=Button(BPSFrame, text='Close', padx=0, pady=0, command = self.myHide)
            btn4.grid(row=10, column=1, sticky=N+W+S+E)

            # Residue selection column
            # Create frame for residue content
            BPSFrame.Res = Frame(BPSFrame)

            # Create Colorize button
            self.colorize.set(1)
            col_cbx = Checkbutton(BPSFrame.Res, text="Colorize residues", command=self.ToggleColor, variable=self.colorize, onvalue=1, offvalue=0)
            col_cbx.grid(row=0,column=0, sticky=N+W)

            # Create group widget
            BPSFrame.gr = Pmw.Group(BPSFrame.Res, tag_text='Residues')
            BPSFrame.gr.grid(row=1, column=0, padx=3, pady=3, sticky=N+W+E)

            # Create residues buttons and checkboxes
            BPSFrame.ccb = []
            for i in range(0,len(self.res_states)):
                var = IntVar()
                var.set(1)
                BPSFrame.ccb.append(Button(BPSFrame.gr.interior(), text='', bg=self.res_states[i][1], activebackground=self.res_states[i][1], relief=RAISED, overrelief=RAISED, height=1, width=2, padx=0, pady=0, command=lambda nr=i: self.changecolor(nr)))
                BPSFrame.ccb[i].grid(row=i,column=0, sticky=N+W)
                rs_cbx = Checkbutton(BPSFrame.gr.interior(), text=self.res_states[i][0], variable=var, onvalue="1", offvalue="0")
                rs_cbx.grid(row=i,column=1, sticky=N+W)
                self.res_states[i][2]=var
            undo1 = Button(BPSFrame.Res, text='Undo all', padx=0, pady=0, command = lambda u="all":self.undo(u))
            undo1.grid(row=2,column=0, sticky=N+W+E)
            undo2 = Button(BPSFrame.Res, text='Undo last', padx=0, pady=0, command = lambda u="last":self.undo(u))
            undo2.grid(row=3,column=0, sticky=N+W+E)
            # Define own colors
            self.create_bp_colors()
            # Make a backup session
            cmd.save('bp_undoall.pse', '(all)')

            # create callback to prevent window kill
            BPSFrame.protocol("WM_DELETE_WINDOW", self.myHide)
            BPSFrame.mainloop()
        else:
            if BPSFrame.state() == "normal":
                self.myHide()
            elif BPSFrame.state() == "withdrawn":
                self.myShow()


    def ToggleColor(self):
        if self.colorize.get()==1:
            for i in range(0,len(self.res_states)):
                BPSFrame.ccb[i].config(state=ACTIVE, bg=self.res_states[i][1], activebackground=self.res_states[i][1], relief=RAISED, overrelief=RAISED)
        else:
            for i in range(0,len(self.res_states)):
                BPSFrame.ccb[i].config(state=DISABLED, bg=self.Default_Color , activebackground=self.Default_Color, relief=RAISED, overrelief=RAISED)

    def changecolor(self,nr):
        new_color = tkColorChooser.askcolor(self.res_states[nr][1])
        rgb_color = map(lambda color:  '%.4f'%float(color/float(65536)) ,new_color[0])
        self.res_states[nr][1]=new_color[1]
        print rgb_color,new_color[1]
        cmd.set_color('bp_plugin_color_'+str(nr),new_color[0],0)
        BPSFrame.ccb[nr].config(bg=new_color[1], activebackground=new_color[1])

    def create_bp_colors(self):
        for i in range(0,len(self.res_states)):
            rgb_color = map(lambda color:  '%.4f'%float(color/float(65536)) ,BPSFrame.winfo_rgb(self.res_states[i][1]))
            cmd.set_color('bp_plugin_color_'+str(i),rgb_color)

    def state(self):
        return map(lambda var: var[2].get(), self.res_states)

    def undo(self,u):
        if u=="all":
            cmd.load('bp_undoall.pse')
        if u=="last":
            cmd.load('bp_undolast.pse')
        self.TogglePlane()

    def scaling(self,v):
        self.TogglePlane()

    def myShow(self):
        self.refresh_list()
        BPSFrame.deiconify()

    def myHide(self):
        cmd.delete("plane_preview")
        BPSFrame.withdraw()

    def get_selections(self):
        BPSurface.selections = []
        for item in cmd.get_names("all"):
            if cmd.get_type(item)=="object:molecule":
                BPSurface.selections.append(item)
            if cmd.get_type(item)=="selection":
                if item[0]<>"_":
                    BPSurface.selections.append(item)

    def refresh_list(self):
        if BPSurface.pp.get()==1:
            self.CreatePlane()
        # Save old values
        old = BPSurface.selections
        try:
            old_cb1=BPSFrame.cb1.getvalue()[0]
        except:
            old_cb1=''
        try:
            old_cb2=BPSFrame.cb2.getvalue()[0]
        except:
            old_cb2=''
        self.get_selections()

        # Update selectionlist if needed
        if not(old==BPSurface.selections):
            BPSFrame.cb1.setlist(BPSurface.selections)
            BPSFrame.cb2.setlist(BPSurface.selections)

        # Is ligand selection still available?
        if not(old_cb1 in BPSurface.selections):
            # No - clear the combobox entry
            BPSFrame.cb1.component('entryfield').clear()
        else:
            # Yes - set the old value
            BPSFrame.cb1.selectitem(old_cb1,setentry=1)

        # Is binding pocket selection still available?
        if not(old_cb2 in BPSurface.selections):
            # No - clear the combobox entry
            BPSFrame.cb2.component('entryfield').clear()
        else:
            # Yes - set the old value
            BPSFrame.cb2.selectitem(old_cb2,setentry=1)

    def clear_flags(self,r,d):
        cmd.flag('ignore','all','clear')
        if r==1:
            cmd.rebuild('all')
        if d==1:
            cmd.delete('bp_cutoff')
            cmd.delete('bp_surface')
            cmd.delete('bp_hydrophobic')
            cmd.delete('bp_aromatic')
            cmd.delete('bp_negative')
            cmd.delete('bp_positive')
            cmd.delete('bp_proline')
            cmd.delete('bp_cysteine')
            cmd.delete('bp_polar')

    def CreatePlane(self):
            view = cmd.get_view()
            # choose the size of the plane
            plane_size = abs(view[11])/4.0
            obj = []

            # Create a plane in camera space
            plane = [
                [ -plane_size,  plane_size, BPSFrame.slider1.get() ],
                [  plane_size,  plane_size, BPSFrame.slider1.get() ],
                [ -plane_size, -plane_size, BPSFrame.slider1.get() ],
                [  plane_size, -plane_size, BPSFrame.slider1.get() ]]

            # Using custom plane?
            if BPSurface.cp.get()==1:
                sinx = sin(pi*BPSFrame.slider2.get()/180)
                siny = sin(pi*BPSFrame.slider3.get()/180)
                cosx = cos(pi*BPSFrame.slider2.get()/180)
                cosy = cos(pi*BPSFrame.slider3.get()/180)
                # Map XY axis rotation
                # cosY  sinX*sinY   -cosX*sinY
                # 0     cosX        sinX
                # sinY  -sinX*cosY  cosX*cosY
                plane = map( lambda p: [
                cosy*p[0] + sinx*siny*p[1] + (-cosx*siny*p[2]),
                cosx * p[1] + sinx* p[2],
                siny * p[0] + (-sinx*cosy*p[1]) + cosx*cosy*p[2]
                ], plane )

            normal = [ 0.0, 0.0, 1.0 ]

            # Transform plane coordinates into model space
            plane = map( lambda p,v=view: [
                v[0] * p[0] + v[1] * p[1] + v[2]* p[2],
                v[3] * p[0] + v[4] * p[1] + v[5]* p[2],
                v[6] * p[0] + v[7] * p[1] + v[8]* p[2]
            ], plane )

            normal = apply( lambda p,v=view:[
                v[0] * p[0] + v[1] * p[1] + v[2]* p[2],
                v[3] * p[0] + v[4] * p[1] + v[5]* p[2],
                v[6] * p[0] + v[7] * p[1] + v[8]* p[2]
            ], (normal,) )

            # Transform position relative to the camera
            plane = map( lambda p,v=view: [
                p[0] + v[9 ] + v[12],
                p[1] + v[10] + v[13],
                p[2] +       + v[14],
            ], plane )
            obj.extend( [ COLOR, 0.5, 0.5, 1.0 ] ) # blueish
            obj.extend( [ BEGIN, TRIANGLE_STRIP ] )
            obj.append( NORMAL )
            obj.extend( normal )
            # draw the plane
            for a in plane:
                obj.append( VERTEX)
                obj.extend(a)
            obj.append( END )
            # delete existing object (if any)
            cmd.delete("plane_preview")
            # now load the new object without zooming
            cmd.set('auto_zoom', 0, quiet=1)
            auto_zoom = cmd.get('auto_zoom')
            cmd.load_cgo(obj,'plane_preview')
            cmd.set('auto_zoom', auto_zoom, quiet=1)

    def TogglePlane(self):
        if BPSurface.pp.get()==0:
            cmd.delete("plane_preview")
        if BPSurface.pp.get()==1:
            self.CreatePlane()

    def ToggleCustomPlane(self):
        self.TogglePlane()
        if BPSurface.cp.get()==0:
            BPSFrame.slider2.grid_remove()
            BPSFrame.slider3.grid_remove()
        if BPSurface.cp.get()==1:
            BPSFrame.slider2.grid(row=3, column=0, columnspan=2, sticky=N+W+S+E,)
            BPSFrame.slider3.grid(row=4, column=0, columnspan=2, sticky=N+W+S+E,)

    def ToggleResidueSelect(self):
        if BPSurface.rv.get()==1:
            BPSFrame.Res.grid(row=0, column=2, rowspan=10, sticky=N+W+E)
        if BPSurface.rv.get()==0:
            BPSFrame.Res.grid_remove()
    
    def ok(self):
        # Has cb1 a value?
        sel1 = BPSFrame.cb1.getvalue()
        if len(sel1) == 0:
            dialog = Pmw.MessageDialog(BPSFrame,
                title = 'Error',
                defaultbutton = 0,
                buttons = ('OK',),
                message_text = 'ERROR: No ligand selection. Please\n choose one from the dropdownlist.'
            )
            dialog.iconname('Selection Error')
            dialog.bell()
            dialog.activate()
            return 0

        # Has cb2 a value?
        sel2 = BPSFrame.cb2.getvalue()
        if len(sel2) == 0:
            dialog = Pmw.MessageDialog(BPSFrame,
                title = 'Error',
                defaultbutton = 0,
                buttons = ('OK',),
                message_text = 'ERROR: No binding pocket selection. Please\n choose one from the dropdownlist.'
            )
            dialog.iconname('Selection Error')
            dialog.bell()
            dialog.activate()
            return 0

        # Are cb1 and cb2 the same?
        if sel1[0]==sel2[0]:
            dialog = Pmw.MessageDialog(BPSFrame,
                title = 'Error',
                defaultbutton = 0,
                buttons = ('OK',),
                message_text = 'ERROR: Ligand and binding pocket must have different selections'
            )
            dialog.iconname('Selection Error')
            dialog.bell()
            dialog.activate()
            return 0

        # Set orgin and normal vector
        cmd.origin(sel1[0])

        vector = [[  0,  0,  1]]

        # Using custom plane?
        if BPSurface.cp.get()==1:
            sinx = sin(pi*BPSFrame.slider2.get()/180)
            siny = sin(pi*BPSFrame.slider3.get()/180)
            cosx = cos(pi*BPSFrame.slider2.get()/180)
            cosy = cos(pi*BPSFrame.slider3.get()/180)
            # Remap the normal vector to custom
            vector = map( lambda p: [
            cosy*p[0] + sinx*siny*p[1] + (-cosx*siny*p[2]),
            cosx * p[1] + sinx* p[2],
            siny * p[0] + (-sinx*cosy*p[1]) + cosx*cosy*p[2]
            ], vector )

        # Remap the normal vector to the model space
        d=cmd.get_view(0)
        vector = map( lambda p,v=d: [
            v[0] * p[0] + v[1] * p[1] + v[2]* p[2],
            v[3] * p[0] + v[4] * p[1] + v[5]* p[2],
            v[6] * p[0] + v[7] * p[1] + v[8]* p[2]
            ], vector )

        # Calculate plane coefficients
        # Ax + By + Cz + D = 0
        # D = -Ax -By -Cz
        A = vector[0][0]
        B = vector[0][1]
        C = vector[0][2]
        D = -vector[0][0]*d[12] - vector[0][1]*d[13] - vector[0][2]*d[14]

        atoms = cmd.get_model(sel2[0])

        # Surface stuff
        self.clear_flags(0,1)
        for at in atoms.atom:
            if (A*(at.coord[0])+B*at.coord[1]+C*at.coord[2]+D-BPSFrame.slider1.get()>0):
                if not('bp_cutoff' in cmd.get_names("selections")):
                    cmd.select('bp_cutoff','index '+str(at.index))
                else:
                    cmd.select('bp_cutoff','bp_cutoff or index '+str(at.index))
        if ('bp_cutoff' in cmd.get_names("selections")):
            if BPSurface.cv.get()==0:
                cmd.flag('ignore','bp_cutoff','reset')
                cmd.delete('indicate')
                cmd.deselect()
                cmd.show('surface',sel2[0])
            else:
                self.clear_flags(0,0)
                cmd.hide('surface',sel2[0])
                cmd.set('auto_zoom','off')
                cmd.create('bp_surface',sel2[0]+' and not bp_cutoff')
                cmd.set('auto_zoom','on')
                cmd.show('surface','bp_surface')
            cmd.rebuild(sel2[0])
        # Selections stuff
        if BPSurface.rv.get()==1:           # make residue selections
            if self.state()[0]==1:          # select hydrophobic
                cmd.select('bp_hydrophobic', 'byres ('+sel2[0]+' and resn ala+gly+ile+leu+met+val)')
            if self.state()[1]==1:          # select aromatic
                cmd.select('bp_aromatic', 'byres ('+sel2[0]+' and resn phe+trp+tyr)')
            if self.state()[2]==1:          # select polar
                cmd.select('bp_polar', 'byres ('+sel2[0]+' and resn his+asn+gln+ser+thr)')
            if self.state()[3]==1:          # select positive
                cmd.select('bp_positive', 'byres ('+sel2[0]+' and resn lys+arg)')
            if self.state()[4]==1:          # select negative
                cmd.select('bp_negative', 'byres ('+sel2[0]+' and resn asp+glu)')
            if self.state()[5]==1:          # select cysteine
                cmd.select('bp_cysteine', 'byres ('+sel2[0]+' and resn cys)')
            if self.state()[6]==1:          # select proline
                cmd.select('bp_proline', 'byres ('+sel2[0]+' and resn pro)')
        cmd.deselect()
        cmd.save('bp_undolast.pse', '(all)')
        if (BPSurface.rv.get()==1 and BPSurface.colorize.get()==1): # colorize?
            if self.state()[0]==1:          # hydrophobic
                cmd.color('bp_plugin_color_0','bp_hydrophobic')
            if self.state()[1]==1:          # select aromatic
                cmd.color('bp_plugin_color_1','bp_aromatic')
            if self.state()[2]==1:          # select polar
                cmd.color('bp_plugin_color_2','bp_polar')
            if self.state()[3]==1:          # select positive
                cmd.color('bp_plugin_color_3','bp_positive')
            if self.state()[4]==1:          # select negative
                cmd.color('bp_plugin_color_4','bp_negative')
            if self.state()[5]==1:          # select cysteine
                cmd.color('bp_plugin_color_5','bp_cysteine')
            if self.state()[6]==1:          # select proline
                cmd.color('bp_plugin_color_6','bp_proline')


The major updates of bp_surface.py version 2.0 are:
  • Ligand and binding site can now be either a selection or an object
  • The orientation of the cutting plane now can be customized
  • Binding site residues can be selected by type: hydrophobic, aromatic, polar, charged positive or negative, proline, cysteine
  • Selected residues can be colorized
  • If you screw up the color settings you can undo the last action or recover the state of the binding site when the script was started the first time (use it with caution)
Here's how the new GUI looks in full view. I like it simple so you only see this window if you check the "Use custom plane" and "Make residues selections". Otherwise parts of the window are collapsed. It's also highly customizable and self-explaining. Lets say I would like to color all polar residues in the binding pocket of 1P93 in red, then I have to check the "Make residues selections" and "Colorize residues" checkboxes, uncheck all types of residues except "polar" and set the color to red by clicking the color button to the left of the "polar" checkbox. In addition I found polar contacts with the ligand and binding pocket residues and hided the generated surface. This looks like this:

Same molecule (Dexamethasone) but this time hydrophobic residues were colored orange. Notice that large parts of the binding site surface are hydrophobic, which isn't a surprise since the bounded ligand is a steroid.

I hope this tool will be helpful ]:)