pymol

.pymol/startup/Draw_Protein_Dimensions.py

@@ -0,0 +1,362 @@
+
+'''
+Calculate and display the dimensions of a protein.
+
+
+This is a first version, please use at your own risk!
+
+REQUIREMENTS
+
+numpy (http://numpy.scipy.org) that should be built into the newers versions of Pymol
+
+(c) Pablo Guardado Calvo
+
+Based on "inertia_tensor.py" (c) 2010 by Mateusz Maciejewski
+
+License: MIT
+
+'''
+
+from __future__ import print_function
+
+__author__  = 'Pablo Guardado Calvo'
+__version__ = '0.1'
+__email__   = 'pablo.guardado (at) gmail.com'
+__date__    = '13/08/2015'
+
+###########################################################################################################################################################
+# USAGE
+#
+# The idea behing this script is to calculate an aproximate minimal bounding box to extract the cell dimensions of a protein. To calculate the minimal bounding
+# is not trivial and usually the Axis Aligned Bounding Box (AABB) does not show up the real dimensions of the protein. This script calculates the inertia tensor
+# of the object, extract the eigenvalues and use them to rotate the molecule (using as rotation matrix the transpose of the eigenvalues matrix). The result is that
+# the molecule is oriented with the inertia axis aligned with the cartesian axis. A new Bounding Box is calculated that is called Inertia Axis Aligned Bounding Box
+#(IABB), whose volume is always lower than AABB volume, and in many cases will correspond with the lowest volume. Of course, maybe it exists another Bounding Box
+# with a lower volume (the minimal Bounding Box).
+#
+# As always with these type of things, you have to use at your own risk. I did not try all the possible combinations, but if you find a bug, do
+# not hesitate to contact me (pablo.guardado (at) gmail.com) or try to modify the code for yourself to correct it.
+#
+# To load the script just type:
+#
+# run path-to-the-script/Draw_Protein_Dimensions.py
+#
+# or if you want something more permanent add the previous line to your .pymolrc file
+#
+# The script works just typing:
+#
+# draw_Protein_Dimensions selection
+#
+# This will draw the cell dimensions of your selection based on a IABB. It also generates the IABB box and the inertia axis, you just need to do "show cgo" to display them.
+#
+# You could also try:
+#
+# draw_BB selection
+#
+# This will draw the AABB and IABB boxes with their cell dimensions and show in the command line their volumes, you can compare both of them.
+############################################################################################################################################################
+
+from pymol import cmd, cgo
+from pymol.cgo import *
+import numpy
+from random import randint
+
+def matriz_inercia(selection):
+	'''
+	DESCRIPTION
+
+	The method calculates the mass center, the inertia tensor and the eigenvalues and eigenvectors
+	for a given selection. Mostly taken from inertia_tensor.py
+	'''
+
+	model = cmd.get_model(selection)
+	totmass = 0.0
+	x,y,z = 0,0,0
+	for a in model.atom:
+		m = a.get_mass()
+		x += a.coord[0]*m
+		y += a.coord[1]*m
+		z += a.coord[2]*m
+		totmass += m
+	global cM
+	cM = numpy.array([x/totmass, y/totmass, z/totmass])
+
+
+	I = []
+	for index in range(9):
+		I.append(0)
+
+	for a in model.atom:
+		temp_x, temp_y, temp_z = a.coord[0], a.coord[1], a.coord[2]
+		temp_x -= x
+		temp_y -= y
+		temp_z -= z
+
+		I[0] += a.get_mass() * (temp_y**2 + temp_z**2)
+		I[1] -= a.get_mass() * temp_x * temp_y
+		I[2] -= a.get_mass() * temp_x * temp_z
+		I[3] -= a.get_mass() * temp_x * temp_y
+		I[4] += a.get_mass() * (temp_x**2 + temp_z**2)
+		I[5] -= a.get_mass() * temp_y * temp_z
+		I[6] -= a.get_mass() * temp_x * temp_z
+		I[7] -= a.get_mass() * temp_y * temp_z
+		I[8] += a.get_mass() * (temp_x**2 + temp_y**2)
+
+	global tensor
+	tensor = numpy.array([(I[0:3]), (I[3:6]), (I[6:9])])
+
+	global autoval, autovect, ord_autoval, ord_autovect
+	autoval, autovect = numpy.linalg.eig(tensor)
+	auto_ord = numpy.argsort(autoval)
+	ord_autoval = autoval[auto_ord]
+	ord_autovect_complete = autovect[:, auto_ord].T
+	ord_autovect = numpy.around(ord_autovect_complete, 3)
+
+	return ord_autoval
+
+
+def draw_inertia_axis(selection):
+	'''
+	DESCRIPTION
+
+	This method draw the inertia axis calculated with the method matriz_inercia.
+	'''
+
+	matriz_inercia(selection)
+	axis1 = ord_autovect[0]
+	x1, y1, z1 = cM[0], cM[1], cM[2]
+	x2, y2, z2 = cM[0]+50*axis1[0], cM[1]+50*axis1[1], cM[2]+50*axis1[2]
+	eje1 = [cgo.CYLINDER, x1, y1, z1, x2, y2, z2, 0.6, 1, 0, 0, 1, 0, 0, 0.0]
+	cmd.load_cgo(eje1, 'Inertia_Axis1')
+	axis2 = ord_autovect[1]
+	x3, y3, z3 = cM[0]+40*axis2[0], cM[1]+40*axis2[1], cM[2]+40*axis2[2]
+	eje1 = [cgo.CYLINDER, x1, y1, z1, x3, y3, z3, 0.6, 1, 0.5, 0, 1, 0.5, 0, 0.0]
+	cmd.load_cgo(eje1, 'Inertia_Axis2')
+	axis4 = ord_autovect[2]
+	x4, y4, z4 = cM[0]+30*axis4[0], cM[1]+30*axis4[1], cM[2]+30*axis4[2]
+	eje1 = [cgo.CYLINDER, x1, y1, z1, x4, y4, z4, 0.6, 1, 1, 0, 1, 1, 0, 0.0]
+	cmd.load_cgo(eje1, 'Inertia_Axis3')
+
+def translacion_cM(selection):
+	'''
+	DESCRIPTION
+
+	Translate the center of mass of the molecule to the origin.
+	'''
+	model = cmd.get_model(selection)
+	totmass = 0.0
+	x,y,z = 0,0,0
+	for a in model.atom:
+		m = a.get_mass()
+		x += a.coord[0]*m
+		y += a.coord[1]*m
+		z += a.coord[2]*m
+		totmass += m
+	cM = numpy.array([x/totmass, y/totmass, z/totmass])
+	trans_array = ([1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, -cM[0], -cM[1], -cM[2], 1])
+	model_trans = cmd.transform_selection(selection, trans_array)
+
+
+def rotacion_orig(selection):
+	'''
+	DESCRIPTION
+
+	Find the proper rotation matrix, i.e. the transpose of the matrix formed by the eigenvectors of the inertia tensor
+	'''
+
+	translacion_cM(selection)
+	matriz_inercia(selection)
+	global transf, transf_array, ord_autovect_array, transf_array_print
+	ord_autovect_array = numpy.array([[ord_autovect[0][0], ord_autovect[0][1], ord_autovect[0][2]],
+	                                  [ord_autovect[1][0], ord_autovect[1][1], ord_autovect[1][2]],
+					  [ord_autovect[2][0], ord_autovect[2][1], ord_autovect[2][2]]])
+	if numpy.linalg.det(ord_autovect_array) == -1:
+		ord_autovect_array = numpy.array([[ord_autovect[2][0], ord_autovect[2][1], ord_autovect[2][2]],
+	                                  	  [ord_autovect[1][0], ord_autovect[1][1], ord_autovect[1][2]],
+	 				  	  [ord_autovect[0][0], ord_autovect[0][1], ord_autovect[0][2]]])
+	transf = numpy.transpose(ord_autovect_array)
+	transf_array = numpy.array([transf[0][0], transf[0][1], transf[0][2], 0,
+	                            transf[1][0], transf[1][1], transf[1][2], 0,
+				    transf[2][0], transf[2][1], transf[2][2], 0,
+				    0, 0, 0, 1])
+
+
+def transformar(selection):
+	'''
+	DESCRIPTION
+
+	Rotate the molecule and draw the inertia axis.
+	'''
+
+	rotacion_orig(selection)
+	model_rot = cmd.transform_selection(selection, transf_array, homogenous=0, transpose=1);
+	draw_inertia_axis(selection)
+
+
+
+def draw_AABB(selection):
+        """
+        DESCRIPTION
+        For a given selection, draw the Axes Aligned bounding box around it without padding. Code taken and modified from DrawBoundingBox.py.
+
+        """
+
+ 	AA_original = selection + "_original"
+	model_orig = cmd.create(AA_original, selection)
+
+        ([min_X, min_Y, min_Z],[max_X, max_Y, max_Z]) = cmd.get_extent(AA_original)
+
+	print("The Axis Aligned Bounding Box (AABB) dimensions are (%.2f, %.2f, %.2f)" % (max_X-min_X, max_Y-min_Y, max_Z-min_Z))
+	print("The Axis Aligned Bounding Box (AABB) volume is %.2f A3" % ((max_X-min_X)*(max_Y-min_Y)*(max_Z-min_Z)))
+
+
+        min_X = min_X
+        min_Y = min_Y
+        min_Z = min_Z
+        max_X = max_X
+        max_Y = max_Y
+        max_Z = max_Z
+
+        boundingBox = [
+                LINEWIDTH, float(2),
+
+                BEGIN, LINES,
+                COLOR, float(1), float(1), float(0),
+
+                VERTEX, min_X, min_Y, min_Z,
+                VERTEX, min_X, min_Y, max_Z,
+                VERTEX, min_X, max_Y, min_Z,
+                VERTEX, min_X, max_Y, max_Z,
+                VERTEX, max_X, min_Y, min_Z,
+                VERTEX, max_X, min_Y, max_Z,
+                VERTEX, max_X, max_Y, min_Z,
+                VERTEX, max_X, max_Y, max_Z,
+                VERTEX, min_X, min_Y, min_Z,
+                VERTEX, max_X, min_Y, min_Z,
+                VERTEX, min_X, max_Y, min_Z,
+                VERTEX, max_X, max_Y, min_Z,
+                VERTEX, min_X, max_Y, max_Z,
+                VERTEX, max_X, max_Y, max_Z,
+                VERTEX, min_X, min_Y, max_Z,
+                VERTEX, max_X, min_Y, max_Z,
+                VERTEX, min_X, min_Y, min_Z,
+                VERTEX, min_X, max_Y, min_Z,
+                VERTEX, max_X, min_Y, min_Z,
+                VERTEX, max_X, max_Y, min_Z,
+                VERTEX, min_X, min_Y, max_Z,
+                VERTEX, min_X, max_Y, max_Z,
+                VERTEX, max_X, min_Y, max_Z,
+                VERTEX, max_X, max_Y, max_Z,
+
+		END
+	]
+
+	p0 = '_0'  + str(randint(0, 100))
+	p1 = '_1' + str(randint(0, 100))
+	p2 = '_2' + str(randint(0, 100))
+	p3 = '_3' + str(randint(0, 100))
+	cmd.pseudoatom (pos=[min_X, min_Y, min_Z], object=p0)
+	cmd.pseudoatom (pos=[min_X, min_Y, max_Z], object=p1)
+	cmd.pseudoatom (pos=[min_X, max_Y, min_Z], object=p2)
+	cmd.pseudoatom (pos=[max_X, min_Y, min_Z], object=p3)
+	cmd.distance(None, p0, p3)
+	cmd.distance(None, p0, p2)
+	cmd.distance(None, p0, p1)
+	cmd.hide("nonbonded")
+
+	boxName = "box_AABB_"  + str(randint(0, 100))
+        cmd.load_cgo(boundingBox,boxName)
+        return boxName
+
+
+def draw_IABB(selection):
+        """
+        DESCRIPTION
+        For a given selection, draw the Inertia Axes Aligned bounding box around it without padding. Code taken and modified from DrawBoundingBox.py.
+
+        """
+
+	transformar(selection)
+
+        ([minX, minY, minZ],[maxX, maxY, maxZ]) = cmd.get_extent(selection)
+
+	print("The Inertia Axis Aligned Bounding Box (IABB) dimensions are (%.2f, %.2f, %.2f)" % (maxX-minX, maxY-minY, maxZ-minZ))
+	print("The Inertia Axis Aligned Bounding Box (IABB) volume is %.2f A3" % ((maxX-minX)*(maxY-minY)*(maxZ-minZ)))
+
+
+        minX = minX
+        minY = minY
+        minZ = minZ
+        maxX = maxX
+        maxY = maxY
+        maxZ = maxZ
+
+        boundingBox = [
+                LINEWIDTH, float(2),
+
+                BEGIN, LINES,
+                COLOR, float(1), float(0), float(0),
+
+                VERTEX, minX, minY, minZ,
+                VERTEX, minX, minY, maxZ,
+                VERTEX, minX, maxY, minZ,
+                VERTEX, minX, maxY, maxZ,
+                VERTEX, maxX, minY, minZ,
+                VERTEX, maxX, minY, maxZ,
+                VERTEX, maxX, maxY, minZ,
+                VERTEX, maxX, maxY, maxZ,
+                VERTEX, minX, minY, minZ,
+                VERTEX, maxX, minY, minZ,
+                VERTEX, minX, maxY, minZ,
+                VERTEX, maxX, maxY, minZ,
+                VERTEX, minX, maxY, maxZ,
+                VERTEX, maxX, maxY, maxZ,
+                VERTEX, minX, minY, maxZ,
+                VERTEX, maxX, minY, maxZ,
+                VERTEX, minX, minY, minZ,
+                VERTEX, minX, maxY, minZ,
+                VERTEX, maxX, minY, minZ,
+                VERTEX, maxX, maxY, minZ,
+                VERTEX, minX, minY, maxZ,
+                VERTEX, minX, maxY, maxZ,
+                VERTEX, maxX, minY, maxZ,
+                VERTEX, maxX, maxY, maxZ,
+
+                END
+        ]
+
+ 	p4 = '_4' + str(randint(0, 100))
+	p5 = '_5' + str(randint(0, 100))
+	p6 = '_6' + str(randint(0, 100))
+	p7 = '_7' + str(randint(0, 100))
+	cmd.pseudoatom (pos=[minX, minY, minZ], object=p4)
+	cmd.pseudoatom (pos=[minX, minY, maxZ], object=p5)
+	cmd.pseudoatom (pos=[minX, maxY, minZ], object=p6)
+	cmd.pseudoatom (pos=[maxX, minY, minZ], object=p7)
+	cmd.distance(None, p4, p7)
+	cmd.distance(None, p4, p6)
+	cmd.distance(None, p4, p5)
+	cmd.hide("nonbonded")
+
+	boxName = "box_IABB_" + str(randint(0, 100))
+        cmd.load_cgo(boundingBox,boxName)
+        return boxName
+
+
+def draw_BB(selection):
+	draw_AABB(selection)
+	draw_IABB(selection)
+
+def draw_Protein_Dimensions(selection):
+	draw_IABB(selection)
+	cmd.hide("cgo")
+
+cmd.extend ("draw_Protein_Dimensions", draw_Protein_Dimensions)
+cmd.extend ("draw_BB", draw_BB)
+
+
+
+
+
+
+

.pymol/startup/axes.py

@@ -0,0 +1,68 @@
+from pymol import cmd
+from chempy import cpv
+
+class PutCenterCallback(object):
+    prev_v = None
+
+    def __init__(self, name, corner=0):
+        self.name = name
+        self.corner = corner
+        self.cb_name = cmd.get_unused_name('_cb')
+
+    def load(self):
+        cmd.load_callback(self, self.cb_name)
+
+    def __call__(self):
+        if self.name not in cmd.get_names('objects'):
+            import threading
+            threading.Thread(None, cmd.delete, args=(self.cb_name,)).start()
+            return
+
+        v = cmd.get_view()
+        if v == self.prev_v:
+            return
+        self.prev_v = v
+
+        t = v[12:15]
+
+        if self.corner:
+            vp = cmd.get_viewport()
+            R_mc = [v[0:3], v[3:6], v[6:9]]
+            off_c = [0.15 * v[11] * vp[0] / vp[1], 0.15 * v[11], 0.0]
+            if self.corner in [2,3]:
+                off_c[0] *= -1
+            if self.corner in [3,4]:
+                off_c[1] *= -1
+            off_m = cpv.transform(R_mc, off_c)
+            t = cpv.add(t, off_m)
+
+        z = -v[11] / 30.0
+        m = [z, 0, 0, t[0] / z, 0, z, 0, t[1] / z, 0, 0, z, t[2] / z, 0, 0, 0, 1]
+        cmd.set_object_ttt(self.name, m, homogenous=1)
+
+def axes(name='axes'):
+    '''
+DESCRIPTION
+
+    Puts coordinate axes to the lower left corner of the viewport.
+    '''
+    from pymol import cgo
+
+    cmd.set('auto_zoom', 0)
+
+    w = 0.06 # cylinder width
+    l = 0.75 # cylinder length
+    h = 0.25 # cone hight
+    d = w * 1.618 # cone base diameter
+
+    obj = [cgo.CYLINDER, 0.0, 0.0, 0.0,   l, 0.0, 0.0, w, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0,
+           cgo.CYLINDER, 0.0, 0.0, 0.0, 0.0,   l, 0.0, w, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0,
+           cgo.CYLINDER, 0.0, 0.0, 0.0, 0.0, 0.0,   l, w, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0,
+           cgo.CONE,   l, 0.0, 0.0, h+l, 0.0, 0.0, d, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0,
+           cgo.CONE, 0.0,   l, 0.0, 0.0, h+l, 0.0, d, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0,
+           cgo.CONE, 0.0, 0.0,   l, 0.0, 0.0, h+l, d, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0]
+
+    PutCenterCallback(name, 1).load()
+    cmd.load_cgo(obj, name)
+
+cmd.extend('axes', axes)

.pymol/startup/center_of_mass.py

@@ -0,0 +1,86 @@
+'''
+See more here: http://www.pymolwiki.org/index.php/center_of_mass
+
+DESCRIPTION
+
+   Places a pseudoatom at the center of mass
+
+   Author: Sean Law
+   Michigan State University
+   slaw (at) msu . edu
+
+SEE ALSO
+
+   pseudoatom, get_com
+'''
+
+from __future__ import print_function
+from pymol import cmd
+
+
+def com(selection, state=None, mass=None, object=None, quiet=1, **kwargs):
+    quiet = int(quiet)
+    if (object == None):
+        try:
+            object = cmd.get_legal_name(selection)
+            object = cmd.get_unused_name(object + "_COM", 0)
+        except AttributeError:
+            object = 'COM'
+    cmd.delete(object)
+
+    if (state != None):
+        x, y, z = get_com(selection, mass=mass, quiet=quiet)
+        if not quiet:
+            print("%f %f %f" % (x, y, z))
+        cmd.pseudoatom(object, pos=[x, y, z], **kwargs)
+        cmd.show("spheres", object)
+    else:
+        for i in range(cmd.count_states()):
+            x, y, z = get_com(selection, mass=mass, state=i + 1, quiet=quiet)
+            if not quiet:
+                print("State %d:%f %f %f" % (i + 1, x, y, z))
+            cmd.pseudoatom(object, pos=[x, y, z], state=i + 1, **kwargs)
+            cmd.show("spheres", 'last ' + object)
+
+cmd.extend("com", com)
+
+
+def get_com(selection, state=1, mass=None, quiet=1):
+    """
+ DESCRIPTION
+
+    Calculates the center of mass
+
+    Author: Sean Law
+    Michigan State University
+    slaw (at) msu . edu
+    """
+    quiet = int(quiet)
+
+    totmass = 0.0
+    if mass != None and not quiet:
+        print("Calculating mass-weighted COM")
+
+    state = int(state)
+    model = cmd.get_model(selection, state)
+    x, y, z = 0, 0, 0
+    for a in model.atom:
+        if (mass != None):
+            m = a.get_mass()
+            x += a.coord[0] * m
+            y += a.coord[1] * m
+            z += a.coord[2] * m
+            totmass += m
+        else:
+            x += a.coord[0]
+            y += a.coord[1]
+            z += a.coord[2]
+
+    if (mass != None):
+        return x / totmass, y / totmass, z / totmass
+    else:
+        return x / len(model.atom), y / len(model.atom), z / len(model.atom)
+
+cmd.extend("get_com", get_com)
+
+# vi:expandtab:sw=3

.pymol/startup/color_h.py

@@ -0,0 +1,128 @@
+# color_h
+# -------
+ 
+# PyMOL command to color protein molecules according to the Eisenberg hydrophobicity scale
+ 
+#
+# Source: http://us.expasy.org/tools/pscale/Hphob.Eisenberg.html
+# Amino acid scale: Normalized consensus hydrophobicity scale
+# Author(s): Eisenberg D., Schwarz E., Komarony M., Wall R.
+# Reference: J. Mol. Biol. 179:125-142 (1984)
+#
+# Amino acid scale values:
+#
+# Ala:  0.620
+# Arg: -2.530
+# Asn: -0.780
+# Asp: -0.900
+# Cys:  0.290
+# Gln: -0.850
+# Glu: -0.740
+# Gly:  0.480
+# His: -0.400
+# Ile:  1.380
+# Leu:  1.060
+# Lys: -1.500
+# Met:  0.640
+# Phe:  1.190
+# Pro:  0.120
+# Ser: -0.180
+# Thr: -0.050
+# Trp:  0.810
+# Tyr:  0.260
+# Val:  1.080
+#
+# Usage:
+# color_h (selection)
+#
+from pymol import cmd
+ 
+def color_h(selection='all'):
+        s = str(selection)
+	print s
+        cmd.set_color('color_ile',[0.996,0.062,0.062])
+        cmd.set_color('color_phe',[0.996,0.109,0.109])
+        cmd.set_color('color_val',[0.992,0.156,0.156])
+        cmd.set_color('color_leu',[0.992,0.207,0.207])
+        cmd.set_color('color_trp',[0.992,0.254,0.254])
+        cmd.set_color('color_met',[0.988,0.301,0.301])
+        cmd.set_color('color_ala',[0.988,0.348,0.348])
+        cmd.set_color('color_gly',[0.984,0.394,0.394])
+        cmd.set_color('color_cys',[0.984,0.445,0.445])
+        cmd.set_color('color_tyr',[0.984,0.492,0.492])
+        cmd.set_color('color_pro',[0.980,0.539,0.539])
+        cmd.set_color('color_thr',[0.980,0.586,0.586])
+        cmd.set_color('color_ser',[0.980,0.637,0.637])
+        cmd.set_color('color_his',[0.977,0.684,0.684])
+        cmd.set_color('color_glu',[0.977,0.730,0.730])
+        cmd.set_color('color_asn',[0.973,0.777,0.777])
+        cmd.set_color('color_gln',[0.973,0.824,0.824])
+        cmd.set_color('color_asp',[0.973,0.875,0.875])
+        cmd.set_color('color_lys',[0.899,0.922,0.922])
+        cmd.set_color('color_arg',[0.899,0.969,0.969])
+        cmd.color("color_ile","("+s+" and resn ile)")
+        cmd.color("color_phe","("+s+" and resn phe)")
+        cmd.color("color_val","("+s+" and resn val)")
+        cmd.color("color_leu","("+s+" and resn leu)")
+        cmd.color("color_trp","("+s+" and resn trp)")
+        cmd.color("color_met","("+s+" and resn met)")
+        cmd.color("color_ala","("+s+" and resn ala)")
+        cmd.color("color_gly","("+s+" and resn gly)")
+        cmd.color("color_cys","("+s+" and resn cys)")
+        cmd.color("color_tyr","("+s+" and resn tyr)")
+        cmd.color("color_pro","("+s+" and resn pro)")
+        cmd.color("color_thr","("+s+" and resn thr)")
+        cmd.color("color_ser","("+s+" and resn ser)")
+        cmd.color("color_his","("+s+" and resn his)")
+        cmd.color("color_glu","("+s+" and resn glu)")
+        cmd.color("color_asn","("+s+" and resn asn)")
+        cmd.color("color_gln","("+s+" and resn gln)")
+        cmd.color("color_asp","("+s+" and resn asp)")
+        cmd.color("color_lys","("+s+" and resn lys)")
+        cmd.color("color_arg","("+s+" and resn arg)")
+cmd.extend('color_h',color_h)
+ 
+def color_h2(selection='all'):
+        s = str(selection)
+	print s
+	cmd.set_color("color_ile2",[0.938,1,0.938])
+	cmd.set_color("color_phe2",[0.891,1,0.891])
+	cmd.set_color("color_val2",[0.844,1,0.844])
+	cmd.set_color("color_leu2",[0.793,1,0.793])
+	cmd.set_color("color_trp2",[0.746,1,0.746])
+	cmd.set_color("color_met2",[0.699,1,0.699])
+	cmd.set_color("color_ala2",[0.652,1,0.652])
+	cmd.set_color("color_gly2",[0.606,1,0.606])
+	cmd.set_color("color_cys2",[0.555,1,0.555])
+	cmd.set_color("color_tyr2",[0.508,1,0.508])
+	cmd.set_color("color_pro2",[0.461,1,0.461])
+	cmd.set_color("color_thr2",[0.414,1,0.414])
+	cmd.set_color("color_ser2",[0.363,1,0.363])
+	cmd.set_color("color_his2",[0.316,1,0.316])
+	cmd.set_color("color_glu2",[0.27,1,0.27])
+	cmd.set_color("color_asn2",[0.223,1,0.223])
+	cmd.set_color("color_gln2",[0.176,1,0.176])
+	cmd.set_color("color_asp2",[0.125,1,0.125])
+	cmd.set_color("color_lys2",[0.078,1,0.078])
+	cmd.set_color("color_arg2",[0.031,1,0.031])
+        cmd.color("color_ile2","("+s+" and resn ile)")
+        cmd.color("color_phe2","("+s+" and resn phe)")
+        cmd.color("color_val2","("+s+" and resn val)")
+        cmd.color("color_leu2","("+s+" and resn leu)")
+        cmd.color("color_trp2","("+s+" and resn trp)")
+        cmd.color("color_met2","("+s+" and resn met)")
+        cmd.color("color_ala2","("+s+" and resn ala)")
+        cmd.color("color_gly2","("+s+" and resn gly)")
+        cmd.color("color_cys2","("+s+" and resn cys)")
+        cmd.color("color_tyr2","("+s+" and resn tyr)")
+        cmd.color("color_pro2","("+s+" and resn pro)")
+        cmd.color("color_thr2","("+s+" and resn thr)")
+        cmd.color("color_ser2","("+s+" and resn ser)")
+        cmd.color("color_his2","("+s+" and resn his)")
+        cmd.color("color_glu2","("+s+" and resn glu)")
+        cmd.color("color_asn2","("+s+" and resn asn)")
+        cmd.color("color_gln2","("+s+" and resn gln)")
+        cmd.color("color_asp2","("+s+" and resn asp)")
+        cmd.color("color_lys2","("+s+" and resn lys)")
+        cmd.color("color_arg2","("+s+" and resn arg)")
+cmd.extend('color_h2',color_h2)

.pymol/startup/dssp_stride.py

@@ -0,0 +1,981 @@
+""" 2011_04_11: Hongbo Zhu
+      PyMOL plugin for running DSSP and display SS using different colors.
+      The DSSP codes for secondary structure are: 
+   
+      H  -  Alpha helix (4-12) 
+      G  -  3-10 helix 
+      I  -  pi helix 
+      E  -  Strand 
+      B  -  Isolated beta-bridge residue 
+      T  -  Turn 
+      S  -  Bend 
+      -  -  None
+      2011_04_24: add stride
+      STRIDE code (taken from stride.doc) is nearly the same as DSSP
+      H         Alpha helix
+      G         3-10 helix
+      I         PI-helix
+      E         Extended conformation
+      B or b    Isolated bridge
+      T         Turn
+      C         Coil (none of the above)
+"""
+
+# Copyright Notice
+# ================
+# 
+# The PyMOL Plugin source code in this file is copyrighted, but you can
+# freely use and copy it as long as you don't change or remove any of
+# the copyright notices.
+# 
+# ----------------------------------------------------------------------
+#               This PyMOL Plugin is Copyright (C) 2011 by 
+#            Hongbo Zhu <hongbo.zhu.cn at googlemail dot com>
+# 
+#                        All Rights Reserved
+# 
+# Permission to use, copy, modify, distribute, and distribute modified
+# versions of this software and its documentation for any purpose and
+# without fee is hereby granted, provided that the above copyright
+# notice appear in all copies and that both the copyright notice and
+# this permission notice appear in supporting documentation, and that
+# the name(s) of the author(s) not be used in advertising or publicity
+# pertaining to distribution of the software without specific, written
+# prior permission.
+# 
+# THE AUTHOR(S) DISCLAIM ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
+# INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.  IN
+# NO EVENT SHALL THE AUTHOR(S) BE LIABLE FOR ANY SPECIAL, INDIRECT OR
+# CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF
+# USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR
+# OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
+# PERFORMANCE OF THIS SOFTWARE.
+# ----------------------------------------------------------------------
+
+# python lib
+import os
+import sys 
+import platform
+if sys.version_info >= (2,4):
+    import subprocess # subprocess is introduced in python 2.4
+import math
+import random
+import tempfile
+import tkSimpleDialog
+import tkMessageBox
+import tkFileDialog
+import tkColorChooser
+
+import Tkinter
+
+# pymol lib
+try: 
+    from pymol import cmd
+    from pymol.cgo import *
+except ImportError:
+    print 'Warning: pymol library cmd not found.'
+    sys.exit(1)
+    
+# external lib    
+try:
+    import Pmw
+except ImportError:
+    print 'Warning: failed to import Pmw. Exit ...'
+    sys.exit(1)
+    
+VERBOSE = True
+
+#################
+## here we go
+#################
+def __init__(self):
+    """ DSSP and Stride plugin for PyMol
+    """
+    self.menuBar.addmenuitem('Plugin', 'command',
+                             'DSSP & Stride', label = 'DSSP & Stride',
+                             command = lambda s=self : DSSPPlugin(s))
+    
+
+#################
+## GUI related
+#################
+class DSSPPlugin:
+  
+    def __init__(self, app):
+        
+        self.parent = app.root
+        self.dialog = Pmw.Dialog(self.parent,
+                                 buttons = ('Run DSSP',
+                                            'Run Stride',
+                                            'Update Color',
+                                            'Update ss',
+                                            'Exit'),
+                                 title = 'DSSP and Stride Plugin for PyMOL',
+                                 command = self.execute)
+        Pmw.setbusycursorattributes(self.dialog.component('hull'))
+
+        # parameters used by DSSP
+        self.pymol_sel     = Tkinter.StringVar()
+        self.dssp_bin      = Tkinter.StringVar()
+        self.stride_bin    = Tkinter.StringVar()
+        self.dssp_rlt_dict = {}
+        self.stride_rlt_dict = {}
+        self.ss_asgn_prog  = None # which program is used to assign ss
+        
+        self.sel_obj_list = [] # there may be more than one seletion or object
+                               # defined by self.pymol_sel
+                               # treat each selection and object separately
+        if 'DSSP_BIN' not in os.environ and 'PYMOL_GIT_MOD' in os.environ:
+            if sys.platform.startswith('linux') and platform.machine() == 'x86_32':
+                initialdir_dssp = os.path.join(os.environ['PYMOL_GIT_MOD'],"DSSP","i86Linux2","dssp-2")
+                os.environ['DSSP_BIN'] = initialdir_dssp
+            elif sys.platform.startswith('linux') and platform.machine() == 'x86_64':
+                initialdir_dssp = os.path.join(os.environ['PYMOL_GIT_MOD'],"DSSP","ia64Linux2","dssp-2")
+                os.environ['DSSP_BIN'] = initialdir_dssp
+            elif sys.platform.startswith('win'):
+                initialdir_dssp = os.path.join(os.environ['PYMOL_GIT_MOD'],"DSSP","win32","dssp.exe")
+                os.environ['DSSP_BIN'] = initialdir_dssp
+            else:
+                pass
+        if 'DSSP_BIN' in os.environ:
+            if VERBOSE: print 'Found DSSP_BIN in environmental variables', os.environ['DSSP_BIN']
+            self.dssp_bin.set(os.environ['DSSP_BIN'])
+        else:
+            if VERBOSE: print 'DSSP_BIN not found in environmental variables.'
+            self.dssp_bin.set('')
+        if 'STRIDE_BIN' not in os.environ and 'PYMOL_GIT_MOD' in os.environ:
+            if sys.platform.startswith('linux') and platform.machine() == 'x86_32':
+                initialdir_stride = os.path.join(os.environ['PYMOL_GIT_MOD'],"Stride","i86Linux2","stride")
+                os.environ['STRIDE_BIN'] = initialdir_stride
+            elif sys.platform.startswith('linux') and platform.machine() == 'x86_64':
+                initialdir_stride = os.path.join(os.environ['PYMOL_GIT_MOD'],"Stride","ia64Linux2","stride")
+                os.environ['STRIDE_BIN'] = initialdir_stride
+            elif sys.platform.startswith('win'):
+                initialdir_stride = os.path.join(os.environ['PYMOL_GIT_MOD'],"Stride","win32","stride.exe")
+                os.environ['STRIDE_BIN'] = initialdir_stride
+            else:
+                pass
+        if 'STRIDE_BIN' in os.environ:
+            if VERBOSE: print 'Found STRIDE_BIN in environmental variables', os.environ['STRIDE_BIN']
+            self.stride_bin.set(os.environ['STRIDE_BIN'])
+        else:
+            if VERBOSE: print 'STRIDE_BIN not found in environmental variables.'
+            self.stride_bin.set('')
+        
+        # DSSP visualization color
+        # - H        Alpha helix (4-12) 
+        # - G        3-10 helix 
+        # - I        pi helix
+        #
+        # - E        Extended strand 
+        # - B        Isolated beta-bridge residue 
+        #
+        # - T        Turn 
+        # - S        Bend 
+        # - -        None
+        # STRIDE does not have S and None, but it has two more
+        # codes: 'b' same as 'B' and 'C' for coil
+        self.DSSP_SSE_list =   ['H','G','I', 'E','B','T','S','-'] # for the record of the order
+        self.STRIDE_SSE_list = ['H','G','I', 'E','B','b','T','C'] # for the record of the order
+        self.SSE_map = {'H':'H',
+                        'G':'H',
+                        'I':'H',
+                        'E':'S',
+                        'B':'S',
+                        'T':'L',
+                        'S':'L',
+                        '-':'L',
+
+                        'b':'S',
+                        'C':'L'
+                        }
+        self.SSE_name = {
+            'H'  :  'Alpha helix', 
+            'G'  :  '3-10 helix', 
+            'I'  :  'Pi helix', 
+                     
+            'E'  :  'Extended strand',
+            'B'  :  'Isolated beta-bridge',
+                     
+            'T'  :  'Turn', 
+            'S'  :  'Bend', 
+            '-'  :  'None',
+
+            'b'  :  'Isolated beta-bridge',
+            'C'  :  'Coil'
+            }
+
+        self.SSE_col_RGB = {
+            'H':(255,  0,  0),
+            'G':(255,  0,128),
+            'I':(255,170,170),
+
+            'E':(255,255,  0),
+            'B':(153,119, 85),
+            
+            'T':(  0,255,255),
+            'S':(153,255,119),
+            '-':(179,179,179),
+
+            'b':(153,119, 85), # same as 'B'
+            'C':(  0,  0,255)
+            }
+        self.SSE_col = {}
+        for sse in self.SSE_col_RGB.keys():
+            self.SSE_col[sse] = '#%s%s%s' % (hex(self.SSE_col_RGB[sse][0])[2:].zfill(2),
+                                             hex(self.SSE_col_RGB[sse][1])[2:].zfill(2),
+                                             hex(self.SSE_col_RGB[sse][2])[2:].zfill(2))
+
+        self.SSE_res_dict = {}
+        self.SSE_sel_dict = {}
+        
+        w = Tkinter.Label(self.dialog.interior(),
+#                          text = '\nDSSP Plugin for PyMOL\nHongbo Zhu, 2011.\n\nColor proteins according to the secondary structure determined by DSSP.',
+                          text = '\nDSSP and Stride Plugin for PyMOL\nby Hongbo Zhu, 2011\n',
+                          background = 'black', foreground = 'green'
+                          )
+        w.pack(expand = 1, fill = 'both', padx = 10, pady = 5)
+
+        # make a few tabs within the dialog
+        self.notebook = Pmw.NoteBook(self.dialog.interior())
+        self.notebook.pack(fill = 'both', expand=1, padx=10, pady=10)
+
+
+        ######################
+        # Tab : Structure Tab
+        ######################
+        page = self.notebook.add('Structure')       
+        self.notebook.tab('Structure').focus_set()
+        group_struc = Tkinter.LabelFrame(page, text = 'Structure')
+        group_struc.pack(fill='both', expand=True, padx=10, pady=5)
+
+        pymol_sel_ent = Pmw.EntryField(group_struc,
+                                       label_text='PyMOL selection/object:',
+                                       labelpos='wn',
+                                       entry_textvariable=self.pymol_sel
+                                       )
+
+        dssp_bin_ent = Pmw.EntryField(group_struc,
+                                      label_text='DSSP binary:', labelpos='wn',
+                                      entry_textvariable=self.dssp_bin,
+                                      entry_width=20)
+        dssp_bin_but = Tkinter.Button(group_struc, text = 'Browse...',
+                                      command = self.getDSSPBin)
+        
+        stride_bin_ent = Pmw.EntryField(group_struc,
+                                      label_text='Stride binary:', labelpos='wn',
+                                      entry_textvariable=self.stride_bin,
+                                      entry_width=20)
+        stride_bin_but = Tkinter.Button(group_struc, text = 'Browse...',
+                                      command = self.getStrideBin)
+
+        # arrange widgets using grid
+        pymol_sel_ent.grid(sticky='we', row=0, column=0,
+                           columnspan=2, padx=5, pady=5)
+        dssp_bin_ent.grid(sticky='we', row=1, column=0, padx=5, pady=1)
+        dssp_bin_but.grid(sticky='we', row=1, column=1, padx=5, pady=1)
+        stride_bin_ent.grid(sticky='we', row=2, column=0, padx=5, pady=1)
+        stride_bin_but.grid(sticky='we', row=2, column=1, padx=5, pady=1)
+        group_struc.columnconfigure(0, weight=9)
+        group_struc.columnconfigure(1, weight=1)
+
+
+        ######################
+        # Tab : Color Tab
+        ######################
+        # H = alpha helix
+        # G = 3-helix (3/10 helix)
+        # I = 5 helix (pi helix)
+
+        # E = extended strand, participates in beta ladder
+        # B = residue in isolated beta-bridge
+
+        # T = hydrogen bonded turn
+        # S = bend
+        # 
+        #  H        Alpha helix (4-12) 
+        #  G        3-10 helix 
+        #  I        pi helix 
+        #  E        Strand 
+        #  B        Isolated beta-bridge residue 
+        #  T        Turn 
+        #  S        Bend 
+        #  -        None 
+        
+        page = self.notebook.add('Color')
+        
+        group_sse_color = Tkinter.LabelFrame(page, text = 'Secondary Structure Element Color')
+        group_sse_color.grid(sticky='eswn', row=0, column=0, padx=10, pady=5)
+
+        # colors for DSSP surface
+        H_col_lab = Tkinter.Label(group_sse_color, text='Alpha helix (H):')
+        self.H_col_but = Tkinter.Button(group_sse_color,
+                                        bg=self.SSE_col['H'],
+                                        activebackground=self.SSE_col['H'],
+                                        command = self.custermizeHColor)
+        G_col_lab = Tkinter.Label(group_sse_color, text='3-10 helix (G):')
+        self.G_col_but = Tkinter.Button(group_sse_color,
+                                        bg=self.SSE_col['G'],
+                                        activebackground=self.SSE_col['G'],
+                                        command = self.custermizeGColor)
+        I_col_lab = Tkinter.Label(group_sse_color, text='PI helix (I):')
+        self.I_col_but = Tkinter.Button(group_sse_color,
+                                        bg=self.SSE_col['I'],
+                                        activebackground=self.SSE_col['G'],
+                                        command = self.custermizeIColor)
+
+
+        E_col_lab = Tkinter.Label(group_sse_color, text='Extended strand (E):')
+        self.E_col_but = Tkinter.Button(group_sse_color,
+                                        bg=self.SSE_col['E'],
+                                        activebackground=self.SSE_col['E'],
+                                        command = self.custermizeEColor)
+        B_col_lab = Tkinter.Label(group_sse_color, text='Isolated beta-bridge (B):')
+        self.B_col_but = Tkinter.Button(group_sse_color,
+                                        bg=self.SSE_col['B'],
+                                        activebackground=self.SSE_col['B'],
+                                        command = self.custermizeBColor)
+
+        T_col_lab = Tkinter.Label(group_sse_color, text='Turn (T):')
+        self.T_col_but = Tkinter.Button(group_sse_color,
+                                        bg=self.SSE_col['T'],
+                                        activebackground=self.SSE_col['T'],
+                                        command = self.custermizeTColor)
+        S_col_lab = Tkinter.Label(group_sse_color, text='Bend (DSSP S):')
+        self.S_col_but = Tkinter.Button(group_sse_color,
+                                        bg=self.SSE_col['S'],
+                                        activebackground=self.SSE_col['S'],
+                                        command = self.custermizeSColor)
+        N_col_lab = Tkinter.Label(group_sse_color, text='None (DSSP NA):')
+        self.N_col_but = Tkinter.Button(group_sse_color,
+                                        bg=self.SSE_col['-'],
+                                        activebackground=self.SSE_col['-'],
+                                        command = self.custermizeNColor)
+         
+
+        b_col_lab = Tkinter.Label(group_sse_color, text='Isolated beta-bridge (Stride b):')
+        self.b_col_but = Tkinter.Button(group_sse_color,
+                                        bg=self.SSE_col['b'],
+                                        activebackground=self.SSE_col['b'],
+                                        command = self.custermizebColor)
+
+        C_col_lab = Tkinter.Label(group_sse_color, text='Coil (Stride C):')
+        self.C_col_but = Tkinter.Button(group_sse_color,
+                                        bg=self.SSE_col['C'],
+                                        activebackground=self.SSE_col['C'],
+                                        command = self.custermizeCColor)
+
+
+
+        H_col_lab.grid(sticky='e', row=0, column=0, padx=5, pady=3)
+        self.H_col_but.grid(sticky='we', row=0, column=1, padx=5, pady=3)
+        G_col_lab.grid(sticky='e', row=1, column=0, padx=5, pady=3)
+        self.G_col_but.grid(sticky='we', row=1, column=1, padx=5, pady=3)
+        I_col_lab.grid(sticky='e', row=2, column=0, padx=5, pady=3)
+        self.I_col_but.grid(sticky='we', row=2, column=1, padx=5, pady=3)
+
+        E_col_lab.grid(sticky='e', row=0, column=2, padx=5, pady=3)
+        self.E_col_but.grid(sticky='we', row=0, column=3, padx=5, pady=3)
+        B_col_lab.grid(sticky='e', row=1, column=2, padx=5, pady=3)
+        self.B_col_but.grid(sticky='we', row=1, column=3, padx=5, pady=3)
+
+        T_col_lab.grid(sticky='e', row=0, column=4, padx=5, pady=3)
+        self.T_col_but.grid(sticky='we', row=0, column=5, padx=5, pady=3)
+        S_col_lab.grid(sticky='e', row=1, column=4, padx=5, pady=3)
+        self.S_col_but.grid(sticky='we', row=1, column=5, padx=5, pady=3)
+        N_col_lab.grid(sticky='e', row=2, column=4, padx=5, pady=3)
+        self.N_col_but.grid(sticky='we', row=2, column=5, padx=5, pady=3)
+
+        b_col_lab.grid(sticky='e', row=2, column=2, padx=5, pady=3)
+        self.b_col_but.grid(sticky='we', row=2, column=3, padx=5, pady=3)
+
+        C_col_lab.grid(sticky='e', row=3, column=4, padx=5, pady=3)
+        self.C_col_but.grid(sticky='we', row=3, column=5, padx=5, pady=3)
+        ######################
+        # Tab : About Tab
+        ######################
+        page = self.notebook.add('About')
+        group_about = Tkinter.LabelFrame(page, text = 'About DSSP and Stride Plugin for PyMOL')
+        group_about.grid(sticky='we', row=0,column=0,padx=5,pady=3)
+        about_plugin = """ Assign and color secondary structures using DSSP or Stride.
+by Hongbo Zhu <hongbo.zhu.cn .at. googlemail.com>
+Please cite this plugin if you use it in a publication.
+Hongbo Zhu. DSSP and Stride plugin for PyMOL, 2011, BIOTEC, TU Dresden.
+"""
+        
+        label_about = Tkinter.Label(group_about,text=about_plugin)
+        label_about.grid(sticky='we', row=0, column=0, padx=5, pady=10)
+        
+        self.notebook.setnaturalsize()    
+
+        return
+
+
+    def getDSSPBin(self):
+        dssp_bin_fname = tkFileDialog.askopenfilename(
+            title='DSSP Binary', initialdir='',
+            filetypes=[('all','*')], parent=self.parent)
+        if dssp_bin_fname: # if nonempty
+            self.dssp_bin.set(dssp_bin_fname)
+        return
+    
+    def getStrideBin(self):
+        stride_bin_fname = tkFileDialog.askopenfilename(
+            title='Stride Binary', initialdir='',
+            filetypes=[('all','*')], parent=self.parent)
+        if stride_bin_fname: # if nonempty
+            self.stride_bin.set(stride_bin_fname)
+        return
+
+    def runDSSPOneObj(self, one_obj_sel):
+        """ Run DSSP on only one object.
+        
+            @param one_obj_sel: the selection/object involving only one object
+            @param type: string
+        """
+        SSE_res = {
+            'H':{}, 'G':{}, 'I':{},
+            'E':{}, 'B':{},
+            'T':{}, 'S':{}, '-':{}
+            }
+        SSE_sel = {
+            'H':None, 'G':None, 'I':None,
+            'E':None, 'B':None,
+            'T':None, 'S':None, '-':None
+            }
+        
+        pdb_fn = None
+        pdb_os_fh, pdb_fn = tempfile.mkstemp(suffix='.pdb') # file os handle, file name
+        os.close(pdb_os_fh)
+        # DSSP 2.0.4 ignores all residues after 1st TER in the same chain
+        v = cmd.get(name='pdb_use_ter_records') 
+        if v: cmd.set(name='pdb_use_ter_records', value=0) # do not insert TER into the pdb 
+        cmd.save(filename=pdb_fn, selection=one_obj_sel)
+        if v: cmd.set(name='pdb_use_ter_records', value=v) # restore old value
+
+        if VERBOSE:
+            print 'Selection %s saved to %s.' % (one_obj_sel, pdb_fn)
+
+        if pdb_fn is None:
+            print 'WARNING: DSSP has no pdb file to work on!'
+            return None
+
+        print 'Running DSSP for %s ...' % (one_obj_sel,)
+        dssp_sse_dict = {}
+        if sys.version_info >= (2,4):
+            dssp_proc = subprocess.Popen([self.dssp_bin.get(), pdb_fn],
+                                         stdout=subprocess.PIPE,
+                                         stderr=subprocess.STDOUT)
+            dssp_stdout, dssp_stderr = dssp_proc.communicate()
+        else: # use os.system + tempfile
+            dssp_tmpout_os_fh, dssp_tmpout_fn = tempfile.mkstemp(suffix='.dssp')
+            os.close(dssp_tmpout_os_fh)
+            dssp_cmd = '%s %s > %s' % (self.dssp_bin.get(), pdb_fn, dssp_tmpout_fn)
+            os.system(dssp_cmd)
+            fh = open(dssp_tmpout_fn)
+            dssp_stdout = ''.join(fh.readlines())
+            fh.close()           
+
+        sse_started = False
+        for line in dssp_stdout.splitlines():
+            if line.startswith('  #  RESIDUE'):
+                sse_started = True
+                continue
+            elif line.startswith(' !!!'):
+                sse_started = False
+                continue
+            elif sse_started:
+                if len(line) < 10 or line[9] == ' ': continue
+                ch,resname = line[11],line[13]
+                residen,sscode = line[5:11].strip(),line[16] # residen = resnum+icode, col 10 is for icode
+                if sscode == ' ': sscode = '-'
+                k = (ch,resname,residen)
+                dssp_sse_dict[k] = sscode
+
+        self.dssp_rlt_dict[one_obj_sel] = dssp_sse_dict
+        print 'Got SSE for %d residues.' % (len(self.dssp_rlt_dict[one_obj_sel]),)
+
+        # group residues according to their SSE, and chain name
+        for k in self.dssp_rlt_dict[one_obj_sel].keys():
+            #res = '/%s//%s/%d%s/' % (sel,k[0],k[1][1], k[1][2].strip()) # sel name, chain ID, res serial num, icode
+            sse = self.dssp_rlt_dict[one_obj_sel][k]
+            chn = k[0]
+            res = k[2]
+            if VERBOSE: print '(%s) and \"%s\"/%s/  sse=%s' % (str(one_obj_sel),
+                                                               chn.strip(), res, sse)
+            SSE_res[sse].setdefault(chn,[]).append(res)
+            
+        self.SSE_res_dict[one_obj_sel] = SSE_res
+
+        for sse in self.DSSP_SSE_list:
+            sse_sel_name = self.selectSSE(one_obj_sel, sse)
+            SSE_sel[sse] = sse_sel_name
+
+        self.SSE_sel_dict[one_obj_sel] = SSE_sel
+
+        print '\nNumber of residues with SSE element:'
+        for sse in self.DSSP_SSE_list:
+            num = sum([len(SSE_res[sse][chn]) for chn in SSE_res[sse]])
+            print '%20s (%s) : %5d' % (self.SSE_name[sse], sse, num)
+        print
+
+        # clean up pdb_fn and dssp_tmpout_fn created by tempfile.mkstemp()
+        if os.path.isfile(pdb_fn):
+            os.remove(pdb_fn)
+        if sys.version_info < (2,4) and os.path.isfile(dssp_tmpout_fn):
+            os.remove(dssp_tmpout_fn)
+
+        return
+
+
+    def runDSSP(self):
+        """
+            @return: whether DSSP has been executed successfully
+            @rtype: boolean
+        """
+        # delete old results
+        self.sel_obj_list  = [] 
+        self.dssp_rlt_dict = {}
+        self.SSE_res_dict  = {}
+        self.SSE_sel_dict  = {}
+        
+        pdb_fn = None
+        sel_name= None
+        sel = self.pymol_sel.get()
+        
+        if len(sel) > 0: # if any pymol selection/object is specified
+            # save the pymol selection/object in the tmp dir
+            all_sel_names = cmd.get_names('all') # get names of all selections            
+            if sel in all_sel_names:
+                if cmd.count_atoms(sel) == 0:
+                    err_msg = 'ERROR: The selection %s is empty.' % (sel,)
+                    print 'ERROR: %s' % (err_msg,)
+                    tkMessageBox.showinfo(title='ERROR', message=err_msg)
+                    return False
+                else:
+                    sel_name = sel
+            # no selection/object with the input name is found
+            # we assume either a single-word selector or
+            # some other selection-expression is used
+            # NOTE: if more than one selection is selected, they should be 
+            #       saved separately and SSE should be calculated for each 
+            #       of the selections.
+            else:
+                print 'The selection/object you specified is not found.'
+                print 'Your input will be interpreted as a selection-expression.'
+                # check whether the selection is empty
+                tmpsel = self.randomSeleName(prefix='your_sele_')
+                cmd.select(tmpsel,sel)
+                if cmd.count_atoms(tmpsel) == 0:
+                    cmd.delete(tmpsel)
+                    err_msg = 'ERROR: The selection %s is empty.' % (sel,)
+                    print 'ERROR: %s' % (err_msg,)
+                    tkMessageBox.showinfo(title='ERROR', message=err_msg)
+                    return False
+                else:
+                    sel_name = tmpsel
+        else:   # what structure do you want DSSP to work on?
+            err_msg = 'No PyMOL selection/object specified!'
+            print 'ERROR: %s' % (err_msg,)
+            tkMessageBox.showinfo(title='ERROR', message=err_msg)
+            return False
+
+        # each object in the selection is treated as an independent struc
+        objlist = cmd.get_object_list(sel_name)
+        self.ss_asgn_prog = 'DSSP'
+        print 'Starting %s ...' % (self.ss_asgn_prog, )
+
+        for objname in objlist:
+            self.sel_obj_list.append('%s and %s' % (sel_name, objname))
+            self.runDSSPOneObj(self.sel_obj_list[-1])
+
+##        cmd.delete(tmpsel)
+
+        return True
+
+
+    def runStrideOneObj(self, one_obj_sel):
+        """ Run Stride on only one object.
+        
+            @param one_obj_sel: the selection/object involving only one object
+            @param type: string
+        """
+        SSE_res = {
+            'H':{}, 'G':{}, 'I':{},
+            'E':{}, 'B':{}, 'b':{},
+            'T':{}, 'C':{}
+            }
+        SSE_sel = {
+            'H':None, 'G':None, 'I':None,
+            'E':None, 'B':None, 'b':None,
+            'T':None, 'C':None
+            }
+        
+        pdb_fn = None
+        pdb_os_fh, pdb_fn = tempfile.mkstemp(suffix='.pdb') # file os handle, file name
+        os.close(pdb_os_fh)
+        cmd.save(filename=pdb_fn, selection=one_obj_sel)
+        if VERBOSE:
+            print 'Selection %s saved to %s.' % (one_obj_sel, pdb_fn)
+
+        if pdb_fn is None:
+            print 'WARNING: Stride has no pdb file to work on!'
+            return None
+        
+        print 'Running Stride for %s ...' % (one_obj_sel,)        
+        stride_sse_dict = {}
+        if sys.version_info >= (2,4):
+            stride_proc = subprocess.Popen([self.stride_bin.get(), pdb_fn],
+                                           stdout=subprocess.PIPE,
+                                           stderr=subprocess.STDOUT
+                                           )
+            stride_stdout, stride_stderr = stride_proc.communicate()
+        else: # use os.system + tempfile
+            stride_tmpout_os_fh, stride_tmpout_fn = tempfile.mkstemp(suffix='.stride')
+            os.close(stride_tmpout_os_fh)
+            stride_cmd = '%s %s > %s' % (self.stride_bin.get(), pdb_fn, stride_tmpout_fn)
+            os.system(stride_cmd)
+            fh = open(stride_tmpout_fn)
+            stride_stdout = ''.join(fh.readlines())
+            fh.close()
+
+        for line in stride_stdout.splitlines():
+            if line.startswith('ASG'):
+                resname,ch=line[5:8], line[9]
+                # according to stride doc, col 12-15 are used for residue number
+                # actually, resnum+icode is used in 12-15.
+                # In addition, if residue number is 4-digit or longer,
+                # the field 12-16 or more are used. 
+
+                if line[15] == ' ': # col 16 is not occupied
+                    residen, sscode = line[11:15].strip(),line[24] # residen = resnum+icode
+                else:
+                    shift = line[15:].find(' ') # find where residen stops
+                    residen, sscode = line[11:15+shift].strip(),line[24+shift]
+
+                k = (ch,resname,residen)
+                stride_sse_dict[k] = sscode
+
+        self.stride_rlt_dict[one_obj_sel] = stride_sse_dict
+        print 'Got SSE for %d residues.' % (len(self.stride_rlt_dict[one_obj_sel]),)
+
+        # group residues according to their SSE, and chain name
+        for k in self.stride_rlt_dict[one_obj_sel].keys():
+            sse=self.stride_rlt_dict[one_obj_sel][k]
+            chn = k[0] # this can be a space!
+            res = k[2]
+            if VERBOSE: print '(%s) and \"%s\"/%s/  sse=%s' % (str(one_obj_sel),
+                                                               chn.strip(), res, sse)
+            SSE_res[sse].setdefault(chn,[]).append(res)
+
+        self.SSE_res_dict[one_obj_sel] = SSE_res
+
+        for sse in self.STRIDE_SSE_list:
+            sse_sel_name = self.selectSSE(one_obj_sel, sse)
+            SSE_sel[sse] = sse_sel_name
+        self.SSE_sel_dict[one_obj_sel] = SSE_sel
+        
+        print '\nNumber of residues with SSE element:'
+        for sse in self.STRIDE_SSE_list:
+            num = sum([len(SSE_res[sse][chn]) for chn in SSE_res[sse]])
+            print '%20s (%s) : %5d' % (self.SSE_name[sse], sse, num)
+        print
+        
+        # clean up pdb_fn and dssp_tmpout_fn created by tempfile.mkstemp()
+        if os.path.isfile(pdb_fn):
+            os.remove(pdb_fn)
+        if sys.version_info < (2,4) and os.path.isfile(stride_tmpout_fn):
+            os.remove(stride_tmpout_fn)        
+        
+        return True
+    
+
+    def runStride(self):
+        """
+        """
+        # delete old results
+        self.sel_obj_list = [] 
+        self.stride_rlt_dict = {}
+        self.SSE_res_dict = {}
+        self.SSE_sel_dict = {}
+
+        pdb_fn = None
+        sel_name= None
+        sel = self.pymol_sel.get()
+        
+        if len(sel) > 0: # if any pymol selection/object is specified
+            all_sel_names = cmd.get_names('all') # get names of all selections            
+            if sel in all_sel_names:
+                if cmd.count_atoms(sel) == 0:
+                    err_msg = 'ERROR: The selection %s is empty.' % (sel,)
+                    print 'ERROR: %s' % (err_msg,)
+                    tkMessageBox.showinfo(title='ERROR', message=err_msg)
+                    return False
+                else:
+                    sel_name = sel
+            # no selection/object with the input name is found
+            # we assume either a single-word selector or
+            # some other selection-expression is uesd
+            else:
+                print 'The selection/object you specified is not found.'
+                print 'Your input will be interpreted as a selection-expression.'
+                tmpsel = self.randomSeleName(prefix='your_sele_')
+                cmd.select(tmpsel,sel)
+                if cmd.count_atoms(tmpsel) == 0:
+                    cmd.delete(tmpsel)
+                    err_msg = 'ERROR: The selection %s is empty.' % (sel,)
+                    print 'ERROR: %s' % (err_msg,)
+                    tkMessageBox.showinfo(title='ERROR', message=err_msg)
+                    return False
+                else:
+                    sel_name = tmpsel
+                
+        else:   # what structure do you want Stride to work on?
+            err_msg = 'No PyMOL selection/object specified!'
+            print 'ERROR: %s' % (err_msg,)
+            tkMessageBox.showinfo(title='ERROR', message=err_msg)
+            return False
+        
+        # each object in the selection is treated as an independent struc
+        objlist = cmd.get_object_list(sel_name)
+        self.ss_asgn_prog = 'Stride'
+        print 'Starting %s ...' % (self.ss_asgn_prog, )
+
+        for objname in objlist:
+            self.sel_obj_list.append('%s and %s' % (sel_name, objname))
+            self.runStrideOneObj(self.sel_obj_list[-1])
+
+        return True
+
+
+    def randomSeleName(self, prefix='sele',suffix=''):
+        """ generate a random selection name.
+        """
+        sel_list = cmd.get_names('all')
+        sel_dict = dict(zip(sel_list, range(len(sel_list))))
+        sel_name = '%s%d%s' % (prefix, random.randint(1000,9999), suffix)
+        while(sel_name in sel_dict):
+            sel_name = '%s%d%s' % (prefix, random.randint(1000,9999), suffix)
+            
+        return sel_name
+        
+
+    def selectSSE(self, sel, sse):
+        """ generate selector for selecting all residues having the given sse.
+            return the selection name.
+        """
+        #sel = self.pymol_sel.get()
+        sel_list_chn = []
+
+        if VERBOSE: print '\nSelecting SSE %s ... \n' % (sse)
+        
+        for chn in self.SSE_res_dict[sel][sse]: # color one chain at a time
+            if chn == ' ': chn_str = '-'
+            else: chn_str = chn
+            if VERBOSE: print 'Selecting SSE %s on chain %s ... \n' % (sse, chn)
+            limit=150  # color every 150 residues
+            sel_name_chn = self.randomSeleName(prefix='%s_%s_%s_' % ('_'.join(sel.split()),chn_str,sse))
+                
+            if len(self.SSE_res_dict[sel][sse][chn]) < limit:
+                #sel_expr = '/%s//%s/%s/' % (sel,chn.strip(), '+'.join(self.SSE_res[sse][chn]))
+                # always quote chain name in case it is empty (otherwise sel misinterpreted)
+                sel_expr = '(%s) and \"%s\"/%s/' % (sel,chn.strip(), '+'.join(self.SSE_res_dict[sel][sse][chn]))
+                cmd.select(sel_name_chn, sel_expr)
+                if VERBOSE: print 'select %s, %s' % (sel_name_chn, sel_expr)
+                sel_list_chn.append(sel_name_chn)
+            else:
+                rn = len(self.SSE_res_dict[sel][sse][chn])
+                print 'total number of res with %s = %d' % (sse,rn)
+                sz = int(math.ceil(rn/float(limit)))
+                sel_list_seg = []
+                for i in xrange(sz):
+                    s,e = i*limit, min((i+1)*limit, rn)
+                    print s,e
+                    #sel_expr = '/%s//%s/%s/' % (sel,chn.strip(), '+'.join(self.SSE_res[sse][chn][s:e]))
+                    # always quote chain name in case it is empty (otherwise sel misinterpreted)
+                    sel_expr = '(%s) and \"%s\"/%s/' % (sel, chn.strip(),
+                                                        '+'.join(self.SSE_res_dict[sel][sse][chn][s:e]))
+                    sel_name_seg =  self.randomSeleName(prefix='%s_%s_%s_tmp_' % ('_'.join(sel.split()),chn_str,sse))
+                    cmd.select(sel_name_seg, sel_expr)
+                    if VERBOSE: print 'select %s, %s' % (sel_name_seg, sel_expr)
+                    sel_list_seg.append(sel_name_seg)
+                    
+                sel_expr = ' or '.join(sel_list_seg)
+                cmd.select(sel_name_chn, sel_expr)
+                if VERBOSE: print 'select %s, %s' % (sel_name_chn, sel_expr)
+                [cmd.delete(asel) for asel in sel_list_seg]
+                sel_list_chn.append(sel_name_chn)
+
+        if len(sel_list_chn) > 0:
+            sel_name = self.randomSeleName(prefix='%s_%s_%s_' % ('_'.join(sel.split()), sse, self.ss_asgn_prog))
+            sel_expr = ' or '.join(sel_list_chn)
+            cmd.select(sel_name,sel_expr)
+            [cmd.delete(asel) for asel in sel_list_chn]
+        else:
+            print 'INFO: No residues are assigned to SSE \'%s\'.' % (sse,)
+            sel_name = None
+        
+        return sel_name
+
+
+    def updateColor(self):
+        if self.ss_asgn_prog is None:
+            err_msg = 'Run DSSP or Stride to assign secondary structures first!'
+            print 'ERROR: %s' % (err_msg,)
+            tkMessageBox.showinfo(title='ERROR', message=err_msg)
+        else:
+            print 'Update color for %s' % (self.pymol_sel.get()),
+            print 'using secondary structure assignment by %s' % (self.ss_asgn_prog,)
+
+            if self.ss_asgn_prog == 'DSSP':     SSE_list = self.DSSP_SSE_list
+            elif self.ss_asgn_prog == 'Stride': SSE_list = self.STRIDE_SSE_list
+
+            for sse in SSE_list: # give color names
+                cmd.set_color('%s_color'%(sse,), self.SSE_col_RGB[sse])
+            for sse in SSE_list: # color each SSE
+                for sel_obj in self.sel_obj_list:
+                    if self.SSE_sel_dict[sel_obj][sse] is not None:
+                        cmd.color('%s_color'%(sse,), self.SSE_sel_dict[sel_obj][sse])
+                        print 'color',self.SSE_sel_dict[sel_obj][sse],',',self.SSE_col_RGB[sse]
+                    else:
+                        print 'No residues with SSE \'%s\' to color.' % (sse,)
+
+        return
+    
+    def updateSS(self):
+        if self.ss_asgn_prog is None:
+            err_msg = 'Run DSSP or Stride to assign secondary structures first!'
+            print 'ERROR: %s' % (err_msg,)
+            tkMessageBox.showinfo(title='ERROR', message=err_msg)
+        else:
+            print 'Update secondary structures for %s' % (self.pymol_sel.get()),
+            print 'using secondary structure assignment by %s' % (self.ss_asgn_prog,)
+            print 'SSE mapping: (H,G,I) ==> H; (E,B,b) ==> S; (T,S,-,C) ==> L'
+
+            if self.ss_asgn_prog == 'DSSP':     SSE_list = self.DSSP_SSE_list
+            elif self.ss_asgn_prog == 'Stride': SSE_list = self.STRIDE_SSE_list
+
+            for sse in SSE_list:
+                for sel_obj in self.sel_obj_list:
+                    if self.SSE_sel_dict[sel_obj][sse] is not None:
+                        cmd.alter(self.SSE_sel_dict[sel_obj][sse], 'ss=\'%s\''% (self.SSE_map[sse],))
+                        print 'alter %s, ss=%s' % (self.SSE_sel_dict[sel_obj][sse], self.SSE_map[sse])
+                    else:
+                        print 'No residue with SSE %s to update ss.' % (sse,)
+
+            # show cartoon for the input selection, and rebuild
+            cmd.show('cartoon',self.pymol_sel.get())
+            cmd.rebuild(self.pymol_sel.get()) 
+            return
+    
+
+    def custermizeHColor(self):
+        self.custermizeSSEColor('H')
+    def custermizeGColor(self):
+        self.custermizeSSEColor('G')
+    def custermizeIColor(self):
+        self.custermizeSSEColor('I')
+
+    def custermizeEColor(self):
+        self.custermizeSSEColor('E')
+    def custermizeBColor(self):
+        self.custermizeSSEColor('B')
+
+    def custermizeTColor(self):
+        self.custermizeSSEColor('T')
+    def custermizeSColor(self):
+        self.custermizeSSEColor('S')
+    def custermizeNColor(self):
+        self.custermizeSSEColor('-')
+
+    def custermizebColor(self):
+        self.custermizeSSEColor('b')
+    def custermizeCColor(self):
+        self.custermizeSSEColor('C')
+
+    def custermizeSSEColor(self,sse):
+        SSE_col_but = {
+            'H':self.H_col_but,
+            'G':self.G_col_but,
+            'I':self.I_col_but,
+
+            'E':self.E_col_but,
+            'B':self.B_col_but,
+            
+            'T':self.T_col_but,
+            'S':self.S_col_but,
+            '-':self.N_col_but,
+            
+            'b':self.b_col_but,
+            'C':self.C_col_but,
+            }
+        try:
+            color_tuple, color = tkColorChooser.askcolor(color=self.SSE_col[sse])
+            if color_tuple is not None and color is not None:
+                self.SSE_col_RGB[sse] = color_tuple
+                self.SSE_col[sse] = color
+                
+                SSE_col_but[sse]['bg']=self.SSE_col[sse]
+                SSE_col_but[sse]['activebackground']=self.SSE_col[sse]
+                SSE_col_but[sse].update()
+        except Tkinter._tkinter.TclError:
+            print 'Old color (%s) will be used.' % (self.mesh_col)
+            
+        
+    def execute(self, butcmd):
+        """ Run the cmd represented by the botton clicked by user.
+        """        
+        if butcmd == 'OK':
+            print 'is everything OK?'
+            
+        elif butcmd == 'Run DSSP':
+            rtn = self.runDSSP()
+            if rtn and VERBOSE: print 'Done with DSSP!'
+            
+        elif butcmd == 'Run Stride':
+            rtn = self.runStride()
+            if rtn and VERBOSE:     print 'Done with Stride!'
+
+        elif butcmd == 'Update Color':
+            self.updateColor()
+            
+        elif butcmd == 'Update ss':
+            self.updateSS()
+
+        elif butcmd == 'Exit':
+            print 'Exiting DSSP and Stride Plugin ...'
+            if __name__ == '__main__':
+                self.parent.destroy()
+            else:
+                self.dialog.withdraw()
+            print 'Done.'
+        else:
+            print 'Exiting DSSP and Stride Plugin because of unknown button click ...'
+            self.dialog.withdraw()
+            print 'Done.'
+            
+    
+    def quit(self):
+        self.dialog.destroy() 
+
+
+#############################################
+#
+#
+# Create demo in root window for testing.
+#
+#
+##############################################
+if __name__ == '__main__':
+    
+    class App:
+        def my_show(self,*args,**kwargs):
+            pass
+
+    app = App()
+    app.root = Tkinter.Tk()
+    Pmw.initialise(app.root)   
+    app.root.title('It seems to work!')
+
+    widget = DSSPPlugin(app)
+    app.root.mainloop()

.pymol/startup/propka.py

@@ -0,0 +1,887 @@
+'''
+Described at PyMOL wiki:
+http://www.pymolwiki.org/index.php/propka
+
+#-------------------------------------------------------------------------------
+# Name:		propka for pymol
+# Purpose:	To fetch and display the pka values for protein of intetest
+#
+# Author:	Troels E. Linnet
+#
+# Created:	14/08/2011
+# Copyright:	(c) Troels E. Linnet 2011
+# Contact:	tlinnet snabela gmail dot com
+# Licence:	Free for all
+#
+#
+#-------------------------------------------------------------------------------
+
+	    The PROPKA method is developed by the
+		  Jensen Research Group
+		 Department of Chemistry
+		University of Copenhagen
+
+	Please cite these references in publications:
+Hui Li, Andrew D. Robertson, and Jan H. Jensen
+"Very Fast Empirical Prediction and Interpretation of Protein pKa Values"
+Proteins, 2005, 61, 704-721.
+
+Delphine C. Bas, David M. Rogers, and Jan H. Jensen
+"Very Fast Prediction and Rationalization of pKa Values for Protein-Ligand Complexes"
+Proteins, 2008, 73, 765-783.
+
+Mats H.M. Olsson, Chresten R. Soendergard, Michal Rostkowski, and Jan H. Jensen
+"PROPKA3: Consistent Treatment of Internal and Surface Residues in Empirical pKa predictions"
+Journal of Chemical Theory and Computation, 2011 7 (2), 525-537
+
+Chresten R. Soendergaard, Mats H.M. Olsson, Michaz Rostkowski, and Jan H. Jensen
+"Improved Treatment of Ligands and Coupling Effects in Empirical Calculation and Rationalization of pKa Values"
+Journal of Chemical Theory and Computation, 2011 in press
+"""
+#-------------------------------------------------------------------------------
+# The script needs mechanize to run.
+# On windows, it is not easy to make additional modules available for pymol. So put in into your working folder.
+#1)The easy manual way:
+#a)Go to: http://wwwsearch.sourceforge.net/mechanize/download.html
+#b)Download mechanize-0.2.5.zip. http://pypi.python.org/packages/source/m/mechanize/mechanize-0.2.5.zip
+#c)Extract to .\mechanize-0.2.5 then move the in-side folder "mechanize" to your folder with propka.py. The rest of .\mechanize-0.2.5 you don't need.
+#You can also see other places where you could put the "mechanize" folder. Write this in pymol to see the paths where pymol is searching for "mechanize"
+# import sys; print(sys.path)
+
+#-------------------------------------------------------------------------------
+"""
+Example for pymol script to start the functions. For example: trypropka.pml
+Execute with pymol or start pymol and: File->Run->trypropka.pml
+##############################################################################################################################################################################################################################
+
+### Point to your directory with your pdb file and where to save the results
+#cd /homes/linnet/Documents/Speciale/5NT-project/Mutant-construct/predict_reactivity/propka
+cd C:/Users/tlinnet/Documents/My Dropbox/Speciale/5NT-project/Mutant-construct/predict_reactivity/propka
+
+### The fastest method is just to write propka. Then the last pymol molecule is assumed and send to server. verbose=yes makes the script gossip mode.
+import propka
+
+fetch 4ins, async=0
+propka
+#fetch 1hp1, async=0
+#propka logtime=_, resi=5-10.20-30, resn=CYS.ATP.TRP, verbose=yes
+
+### Fetch 4ins from web. async make sure, we dont execute script before molecule is loaded. The resi and resn prints the interesting results right to command line.
+#fetch 4ins, async=0
+#propka chain=*, resi=5-10.20-30, resn=ASP.CYS, logtime=_
+
+### If there is no web connection, one can process a local .pka file. Either from a previous run or from a downloaded propka webpage result.
+### Then run and point to .pka file with: pkafile=./Results_propka/pkafile.pka Remember the dot "." in the start, to make it start in the current directory.
+#load 4ins.pdb
+#propka pkafile=./Results_propka/4ins_.pka, resi=18.25-30, resn=cys,
+
+### Some more examples. This molecule has 550 residues, so takes a longer time. We select to run the last molecule, by writing: molecule=1hp1
+#fetch 4ins, async=0
+#fetch 1hp1, async=0
+#propka molecule=1hp1, chain=A, resi=300-308.513, resn=CYS.ATP.TRP, logtime=_, verbose=no, showresult=no
+#propka molecule=1hp1, pkafile=./Results_propka/1hp1_.pka, verbose=yes
+
+### One can also just make a lookup for a protein. Use function: getpropka
+### Note. This does only print the result to the pymol command line
+#getpropka source=ID, PDBID=4ake, logtime=_, showresult=yes
+#getpropka source=ID, PDBID=4ins, logtime=_, server_wait=10.0, verbose=yes, showresult=no
+############################################Input parameters: propka############################################
+############# The order of input and changable things:
+############# propka(molecule="NIL",chain="*",resi="0",resn="NIL",method="upload",logtime=time.strftime("%m%d",time.localtime()),server_wait=3.0,version="v3.1",verbose="no",showresult="no",pkafile="NIL")
+# method : method=upload is default. This sends .pdb file and request result from propka server.
+## method=file will only process a manual .pka file, and write a pymol command file. No use of mechanize.
+## If one points to an local .pka file, then method is auto-changed to method=file. This is handsome in off-line environment, ex. teaching or seminar.
+# pkafile: Write the path to .pka file. Ex: pkafile=./Results_propka/4ins_.pka
+# molecule : name of the molecule. Ending of file is assumed to be .pdb
+# chain : which chains are saved to file, before molecule file is send to server. Separate with "." Ex: chain=A.b
+# resi : Select by residue number, which residues should be printed to screen and saved to the log file: /Results_propka/_Results.log.
+## Separate with "." or make ranges with "-". Ex: resi=35.40-50
+# resn : Select by residue name, which residues should be printed to screen and saved to the log file: /Results_propka/_Results.log.
+## Separate with "." Ex: resn=cys.tyr
+# logtime : Each execution give a set of files with the job id=logtime. If logtime is not provided, the current time is used.
+## Normal it usefull to set it empty. Ex: logtime=_
+# verbose : Verbose is switch, to turn on messages for the mechanize section. This is handsome to see how mechanize works, and for error searching.
+# showresult : Switch, to turn on all results in pymol command window. Ex: showresult=yes
+# server_wait=10.0 is default. This defines how long time between asking the server for a result. Set no lower than 3 seconds.
+# version=v3.1 is default. This is what version of propka which would be used.
+## Possible: 'v3.1','v3.0','v2.0'. If a newer version is available than the current v3.1, a error message is raised to make user update the script.
+############################################Input parameters: getpropka############################################
+############# The order of input and changable things:
+############# getpropka(PDB="NIL",chain="*",resi="0",resn="NIL",source="upload",PDBID="",logtime=time.strftime("%Y%m%d%H%M%S",time.localtime()),server_wait=3.0,version="v3.1",verbose="no",showresult="no")
+# PDB: points the path to a .pdb file. This is auto-set from propka function.
+# source : source=upload is default and is set at the propka webpage.
+# source=ID, PDBID=4ake , one can print to the command line, the pka value for any official pdb ID. No files are displayed in pymol.
+# PDBID: is used as the 4 number/letter pdb code, when invoking source=ID.
+
+##############################################################################################################################################################################################################################
+'''
+
+from __future__ import print_function
+
+try:
+    from pymol import cmd
+    runningpymol = 'yes'
+except:
+    runningpymol = 'no'
+    pass
+import time
+import platform
+import os
+import sys
+
+
+def propka(molecule="NIL", chain="*", resi="0", resn="NIL", method="upload", logtime=time.strftime("%m%d", time.localtime()), server_wait=3.0, version="v3.1", verbose="no", showresult="no", pkafile="NIL", makebonds="yes"):
+    Script_Version = "20110823"
+    # First we have to be sure, we give reasonable arguments
+    if pkafile != "NIL":
+        method = 'file'
+    assert method in ['upload', 'file'], "'method' has to be either: method=upload or method=file"
+    # If molecule="all", then try to get the last molecule
+    ##assert molecule not in ['NIL'], "You always have to provide molecule name. Example: molecule=4ins"
+    if molecule == "NIL":
+        assert len(cmd.get_names()) != 0, "Did you forget to load a molecule? There are no objects in pymol."
+        molecule = cmd.get_names()[-1]
+    # To print out to screen for selected residues. Can be separated with "." or make ranges with "-". Example: resi="4-8.10"
+    if resi != "0":
+        resi_range = ResiRange(resi)
+    else:
+        resi_range = []
+    # Also works for residue names. They are all converted to bigger letters. Example: resn="cys.Tyr"
+    if resn != "NIL":
+        resn_range = ResnRange(resn)
+    else:
+        resn_range = resn
+    # Make chain range, and upper case.
+    chain = ChainRange(chain)
+    # Make result directory. We also the absolut path to the new directory.
+    Newdir = createdirs()
+    if method == "upload":
+        # We try to load mechanize. If this fail, one can always get the .pka file manual and the run: method=file
+        try:
+            from modules import mechanize
+            importedmechanize = 'yes'
+        except ImportError:
+            print("Import error. Is a module missing?")
+            print(sys.exc_info())
+            print("Look if missing module is in your python path\n%s" % sys.path)
+            importedmechanize = 'no'
+            import modules.mechanize as mechanize
+        # The name for the new molecule
+        newmolecule = "%s%s" % (molecule, logtime)
+        # Create the new molecule from original loaded and for the specified chains. Save it, and disable the old molecule.
+        cmd.create("%s" % newmolecule, "%s and chain %s" % (molecule, chain))
+        cmd.save("%s%s.pdb" % (Newdir, newmolecule), "%s" % newmolecule)
+        cmd.disable("%s" % molecule)
+        if molecule == "all":
+            cmd.enable("%s" % molecule)
+            cmd.show("cartoon", "%s" % molecule)
+        # Let the new molecule be shown in cartoon.
+        cmd.hide("everything", "%s" % newmolecule)
+        cmd.show("cartoon", "%s" % newmolecule)
+        # Make the absolut path to the newly created .pdb file.
+        PDB = "%s%s.pdb" % (Newdir, newmolecule)
+        source = "upload"
+        PDBID = ""
+        # Request server, and get the absolut path to the result file.
+        pkafile = getpropka(PDB, chain, resi, resn, source, PDBID, logtime, server_wait, version, verbose, showresult)
+        # Open the result file and put in into a handy list.
+        list_results, ligands_results = importpropkaresult(pkafile)
+    if method == "file":
+        assert pkafile not in ['NIL'], "You have to provide path to file. Example: pkafile=./Results_propka/4ins_2011.pka"
+        assert ".pka" in pkafile, 'The propka result file should end with ".pka" \nExample: pkafile=./Results_propka/4ins_2011.pka \npkafile=%s' % (pkafile)
+        # The name for the molecule we pass to the writing script of pymol commands
+        newmolecule = "%s" % molecule
+        cmd.hide("everything", "%s" % newmolecule)
+        cmd.show("cartoon", "%s" % newmolecule)
+        # We open the result file we have got in the manual way and put in into a handy list.
+        list_results, ligands_results = importpropkaresult(pkafile)
+        # Then we print the interesting residues to the screen.
+        printpropkaresult(list_results, resi, resi_range, resn, resn_range, showresult, ligands_results)
+    # Now create the pymol command file. This should label the protein. We get back the absolut path to the file, so we can execute it.
+    result_pka_pymol_name = writepymolcmd(newmolecule, pkafile, verbose, makebonds)
+    # Now run our command file. But only if we are running pymol.
+    if runningpymol == 'yes':
+        cmd.do("run %s" % result_pka_pymol_name)
+    ##if runningpymol=='yes': cmd.do("@%s"%result_pka_pymol_name)
+    return(list_results)
+if runningpymol != 'no':
+    cmd.extend("propka", propka)
+
+
+def getpropka(PDB="NIL", chain="*", resi="0", resn="NIL", source="upload", PDBID="", logtime=time.strftime("%Y%m%d%H%M%S", time.localtime()), server_wait=3.0, version="v3.1", verbose="no", showresult="no"):
+    try:
+        import modules.mechanize as mechanize
+        importedmechanize = 'yes'
+    except ImportError:
+        print("Import error. Is a module missing?")
+        print(sys.exc_info())
+        print("Look if missing module is in your python path \n %s" % sys.path)
+        importedmechanize = 'no'
+    propka_v_201108 = 3.1
+    url = "http://propka.ki.ku.dk/"
+    assert version in ['v2.0', 'v3.0', 'v3.1'], "'version' has to be either: 'v2.0', 'v3.0', 'v3.1'"
+    assert source in ['ID', 'upload', 'addr', 'input_file'], "'source' has to be either: 'ID', 'upload', 'addr', 'input_file'"
+    if source == "upload":
+        assert PDB not in ['NIL'], "You always have to provide PDB path. Example: PDB=.\Results_propka\4ins2011.pdb"
+    if source == "ID":
+        assert len(PDBID) == 4, "PDBID has to be 4 characters"
+    # To print out to screen for selected residues. Can be separated with "." or make ranges with "-". Example: resi="4-8.10"
+    if resi != "0":
+        resi_range = ResiRange(resi)
+    else:
+        resi_range = []
+    # Also works for residue names. They are all converted to bigger letters. Example: resn="cys.Tyr"
+    if resn != "NIL":
+        resn_range = ResnRange(resn)
+    else:
+        resn_range = resn
+    # Start the browser
+    br = mechanize.Browser()
+    # We pass to the server, that we are not a browser, but this python script. Can be used for statistics at the propka server.
+    br.addheaders = [('User-agent', 'pythonMechanizeClient')]
+    # To turn on debugging messages
+    # br.set_debug_http(True)
+    # To open the start page.
+    page_start = br.open(url)
+    read_start = page_start.read()
+    if verbose == 'yes':
+        print(br.title())
+        print(br.geturl())
+    # To get available forms
+    page_forms = [f.name for f in br.forms()]
+    if verbose == 'yes':
+        print(page_forms)
+    # Select first form
+    br.select_form(name=page_forms[0])
+    # Print the current selected form, so we see that we values we start with.
+    if verbose == 'yes':
+        print(br.form)
+    # Print the parameters of the 'version' RadioControl button and current value
+    if verbose == 'yes':
+        print(br.find_control(name='version'), br.find_control(name='version').value)
+    # This is to check, that the current script is "up-to-date".
+    propka_v_present = float(br.find_control(name='version').value[0].replace('v', ''))
+    if propka_v_present > propka_v_201108:
+        raise UserWarning('\nNew version of propka exist.\nCheck/Update your script.\nPresent:v%s > Script:v%s' % (propka_v_present, propka_v_201108))
+    # Change the parameters of the 'version' radio button and then reprint the new value. Input has to be in a list [input].
+    br.form['version'] = [version]
+    if verbose == 'yes':
+        print(br.find_control(name='version').value)
+    # Print the parameters of the 'source' RadioControl button and current value
+    if verbose == 'yes':
+        print(br.find_control(name='source'), br.find_control(name='source').value)
+    # Change the parameters of the 'source' radio button and then reprint the new value. Input has to be in a list [input].
+    br.form['source'] = [source]
+    if verbose == 'yes':
+        print(br.find_control(name='source').value)
+    # This step was the must strange and took a long time. For finding the information and the double negative way.
+    # One have to enable the pdb button. Read more here: http://wwwsearch.sourceforge.net/old/ClientForm/ ("# All Controls may be disabled.....)
+    PDBID_control = br.find_control("PDBID")
+    PDB_control = br.find_control("PDB")
+    if verbose == 'yes':
+        print(PDBID_control.disabled, PDB_control.disabled)
+    if source == "ID":
+        PDBID_control.disabled = False
+        PDB_control.disabled = True
+    if source == "upload":
+        PDBID_control.disabled = True
+        PDB_control.disabled = False
+    if verbose == 'yes':
+        print(PDBID_control.disabled, PDB_control.disabled)
+    # We create the result dir, and take with us the 'path' to the result dir.
+    Newdir = createdirs()
+    # Open all the files, and assign them.
+    if source == "upload":
+        filename = PDB
+    if source == "ID":
+        filename = PDBID
+    files = openfiles(Newdir, filename, logtime, source)
+    result_pka_file = files[0]
+    result_input_pka_file = files[1]
+    result_log = files[2]
+    filepath = files[3]
+    result_pka_pkafile = files[4]
+    result_pka_file_stripped = files[5]
+    result_pka_file_bonds = files[6]
+    # Print the parameters of the 'PDBID' TextControl button and current value
+    if source == "ID" and verbose == 'yes':
+        print(br.find_control(name='PDBID'))
+        print(br.find_control(name='PDBID').value)
+    # Change the parameters of the 'PDBID' TextControl and then reprint the new value. Input has just to be a string.
+    if source == "ID":
+        br.form["PDBID"] = PDBID
+    if source == "ID" and verbose == 'yes':
+        print(br.find_control(name='PDBID').value)
+    # Print the parameters of the 'PDB' TextControl button and current value
+    if source == "upload" and verbose == 'yes':
+        print(br.find_control(name='PDB'))
+        print(br.find_control(name='PDB').value)
+    # Change the parameters of the 'PDB' FileControl and then reprint the new value. Input has just to be a string.
+    if source == "upload":
+        PDBfilename = PDB
+        PDBfilenamepath = PDB
+    if source == "upload":
+        br.form.add_file(open(PDBfilename), 'text/plain', PDBfilenamepath, name='PDB')
+    if source == "upload" and verbose == 'yes':
+        print(br.find_control(name='PDB'))
+        print(br.find_control(name='PDB').value)
+    # Now reprint the current selected form, so we see that we have the right values.
+    if verbose == 'yes':
+        print(br.form)
+    # Make "how" we would like the next request. We would like to "Click the submit button", but we have not opened the request yet.
+    req = br.click(type="submit", nr=0)
+    # Have to pass by a mechanize exception. Thats the reason for the why True
+    # The error was due to: br.open(req)
+    # mechanize._response.httperror_seek_wrapper: HTTP Error refresh: The HTTP server returned a redirect error that would lead to an infinite loop.
+    # The last 30x error message was:
+    # OK
+    # I haven't been able to find the refresh problem or extend the time. So we make a pass on the raised exception.
+    try:
+        print("Now sending request to server")
+        br.open(req)
+    # If there is raised an exception, we jump through to the result page after some sleep.
+    except mechanize.HTTPError:
+        # We can extract the jobid from the current browser url.
+        jobid = br.geturl()[32:-5]
+        # We notice how the script at the server presents the final result page.
+        url_result = url + "pka/" + jobid + ".html"
+        # Now we continue to try to find the result page, until we have succes. If page doesn't exist, we wait a little.
+        while True:
+            print("Result still not there. Waiting %s seconds more" % server_wait)
+            time.sleep(float(server_wait))
+            # To pass the "break" after the exception, we make a hack, wait and then go to the result page, which is the jobid.
+            try:
+                page_result = br.open(url_result)
+                read_result = page_result.read()
+                # If we don't receive a error in getting the result page, we break out of the while loop.
+                break
+            # If the page doesn't exist yet. We go back in the while loop.
+            except mechanize.HTTPError:
+                # Wait another round
+                pass
+            # If we get a timeout, we also wait.
+            except mechanize.URLError:
+                # Wait another round
+                pass
+    htmlresult = "The detailed result is now available at: %s" % br.geturl()
+    print(htmlresult)
+    read_result = br.response().read()
+    # Now save the available links from the current page. But only links that satisfy the expression.
+    links_result = []
+    for l in br.links(url_regex='http://propka.ki.ku.dk/pka'):
+        links_result.append(l)
+    # We also extract the information for neighbour bons. This is given in the url links.
+    bonds = []
+    for l in br.links(url_regex='http://propka.ki.ku.dk/view/new_view.cgi'):
+        l_split = str(l).split()
+        lresn = l_split[2]
+        lresi = l_split[3]
+        lchain = l_split[4]
+        lurl = l_split[1]
+        lurl_split = lurl.split("&")
+        lresn2 = lurl_split[1]
+        lchain2 = lurl_split[2]
+        lpka = lurl_split[3]
+        ldesolvation = lurl_split[4]
+        lneighbours = lurl_split[5:]
+        for i in range(len(lneighbours)):
+            bonds.append([lresn, lresi, lchain, lresn2, lchain2, lpka, ldesolvation, lneighbours[i]])
+    # Now follow the link to the .propka_input resultpage
+    if len(links_result) > 1:
+        br.follow_link(links_result[1])
+    # Now get the page text for the current link
+    if len(links_result) > 1:
+        read_result1 = br.response().read()
+    # Save the result
+    if len(links_result) > 1:
+        result_input_pka_file.write(read_result1)
+    # Now follow the link to the .pka resultpage
+    if len(links_result) > 1:
+        br.back()
+    result_input_pka_file.close()
+    # Now follow first link. "Should be" available for all versions of propka.
+    br.follow_link(links_result[0])
+    # Now get the page for the current link
+    read_result0 = br.response().read()
+    # Save the result and close file.
+    result_pka_file.write(read_result0)
+    result_pka_file.close()
+    # Now get the result in a list, which is sorted
+    list_results, ligands_results = importpropkaresult(result_pka_pkafile)
+    # Print to log file
+    result_log.write("# executed: %s \n# logtime: %s \n# source=%s \n# PDB=%s \n# chain=%s \n# PDBID=%s \n# server_wait=%s version=%s verbose=%s showresult=%s \n# resi=%s resn=%s\n# %s \n" % (time.strftime("%Y%m%d%H%M%S", time.localtime()), logtime, source, PDB, chain, PDBID, server_wait, version, verbose, showresult, resi, resn, htmlresult))
+    # Print to screen
+    printpropkaresult(list_results, resi, resi_range, resn, resn_range, showresult, ligands_results)
+    # Now write to log and the stripped file
+    for l in list_results:
+        if resi != "0" and int(l[1]) in resi_range:
+            result_log.write("%3s %3s %s %6s %3s %5s %3s %4s %s" % (l[0], l[1], l[2], l[3], l[4], l[5], l[6], l[7], l[8]) + '\n')
+        if resn != "NIL" and l[0] in resn_range and int(l[1]) not in resi_range:
+            result_log.write("%3s %3s %s %6s %3s %5s %3s %4s %s" % (l[0], l[1], l[2], l[3], l[4], l[5], l[6], l[7], l[8]) + '\n')
+        result_pka_file_stripped.write("%3s %3s %s %6s %3s %5s %3s %4s %s" % (l[0], l[1], l[2], l[3], l[4], l[5], l[6], l[7], l[8]) + '\n')
+    for l in ligands_results:
+        if resn != "NIL" and l[0] in resn_range:
+            result_log.write("%3s %3s %s %6s %3s %5s %3s %4s %s" % (l[0], l[1], l[2], l[3], l[4], l[5], l[6], l[7], l[8]) + '\n')
+        result_pka_file_stripped.write("%3s %3s %s %6s %3s %5s %3s %4s %s" % (l[0], l[1], l[2], l[3], l[4], l[5], l[6], l[7], l[8]) + '\n')
+    result_pka_file_stripped.close()
+    result_log.close()
+    # Now handle the bonds. We have to delete dublicates first.
+    bonds.sort()
+    last = bonds[-1]
+    for i in range(len(bonds) - 2, -1, -1):
+        if last == bonds[i]:
+            del bonds[i]
+        else:
+            last = bonds[i]
+    # Now make a selection for known residue
+    bonds_selected = []
+    bonds_ligands = []
+    for l in bonds:
+        if l[0][6:] in ['ASP', 'GLU', 'ARG', 'LYS', 'HIS', 'CYS', 'TYR', 'C-', 'N+']:
+            bonds_selected.append(l)
+        else:
+            bonds_ligands.append(l)
+    # And now sort it.
+    bonds_selected.sort(key=lambda residue: int(residue[1]))
+    # Now write it to file
+    bonddic = {'=': ' ', ':': ' ', ',': ' ', "'": " "}
+    for l in bonds_selected:
+        nb = replace_all(l[7], bonddic)
+        result_pka_file_bonds.write("%3s %3s %s %7s %7s %9s %17s %s" % (l[0][6:], l[1], l[2][:1], l[3][8:], l[4], l[5], l[6], nb) + '\n')
+    for l in bonds_ligands:
+        nb = replace_all(l[7], bonddic)
+        result_pka_file_bonds.write("%3s %3s %s %7s %7s %9s %17s %s" % (l[0][6:], l[1], l[2][:1], l[3][8:], l[4], l[5], l[6], nb) + '\n')
+    result_pka_file_bonds.close()
+    return(result_pka_pkafile)
+if runningpymol != 'no':
+    cmd.extend("getpropka", getpropka)
+
+
+def openpymolfiles(pkafile):
+    result_pka_pymol_name = pkafile.replace(".pka", ".pml")
+    result_pka_pymol = open(result_pka_pymol_name, "w")
+    return(result_pka_pymol, result_pka_pymol_name)
+
+
+def printpropkaresult(list_results, resi, resi_range, resn, resn_range, showresult, ligands_results):
+    for l in list_results:
+        if resi != "0" and int(l[1]) in resi_range:
+            if showresult != 'yes':
+                print("%3s %3s %s %6s %3s %5s %3s %4s %s" % (l[0], l[1], l[2], l[3], l[4], l[5], l[6], l[7], l[8]))
+        if resn != "NIL" and l[0] in resn_range and int(l[1]) not in resi_range:
+            if showresult != 'yes':
+                print("%3s %3s %s %6s %3s %5s %3s %4s %s" % (l[0], l[1], l[2], l[3], l[4], l[5], l[6], l[7], l[8]))
+        if showresult == 'yes':
+            print("%3s %3s %s %6s %3s %5s %3s %4s %s" % (l[0], l[1], l[2], l[3], l[4], l[5], l[6], l[7], l[8]))
+    for l in ligands_results:
+        if resn != "NIL" and l[0] in resn_range:
+            if showresult != 'yes':
+                print("%3s %3s %s %6s %3s %5s %3s %4s %s" % (l[0], l[1], l[2], l[3], l[4], l[5], l[6], l[7], l[8]))
+        if showresult == 'yes':
+            print("%3s %3s %s %6s %3s %5s %3s %4s %s" % (l[0], l[1], l[2], l[3], l[4], l[5], l[6], l[7], l[8]))
+
+
+def importpropkaresult(result_pka_pkafile):
+    result_pka_file = open(result_pka_pkafile, "r")
+    list_results = []
+    ligands_results = []
+    ##bonding_partners = []
+    for l in result_pka_file:
+        if not l.strip():
+            continue
+        else:
+            # To search for the right lines
+            if l.strip().split()[0] in ['ASP', 'GLU', 'ARG', 'LYS', 'HIS', 'CYS', 'TYR', 'C-', 'N+'] and len(l.strip().split()) > 20:
+                list_results.append([l.strip().split()[0], l.strip().split()[1], l.strip().split()[2], l.strip().split()[3], l.strip().split()[4], l.strip().split()[6], l.strip().split()[7], l.strip().split()[8], l.strip().split()[9]])
+                # bonding_partners.append(l.strip().split()[11]);bonding_partners.append(l.strip().split()[15]);bonding_partners.append(l.strip().split()[19])
+            if l.strip().split()[0] not in ['ASP', 'GLU', 'ARG', 'LYS', 'HIS', 'CYS', 'TYR', 'C-', 'N+'] and len(l.strip().split()) > 20:
+                ligands_results.append([l.strip().split()[0], l.strip().split()[1], l.strip().split()[2], l.strip().split()[3], l.strip().split()[4], l.strip().split()[6], l.strip().split()[7], l.strip().split()[8], l.strip().split()[9]])
+                # bonding_partners.append(l.strip().split()[11]);bonding_partners.append(l.strip().split()[15]);bonding_partners.append(l.strip().split()[19])
+    # Sort the result after the residue number and then chain.
+    list_results.sort(key=lambda residue: int(residue[1]))
+    list_results.sort(key=lambda chain: chain[2])
+    # bonding_partners=uniqifi(bonding_partners)
+    ##bonding_partners[:] = [x for x in bonding_partners if x != "XXX"]
+    result_pka_file.close()
+    return(list_results, ligands_results)
+
+
+def importpropkabonds(result_pka_pkafile):
+    bonds = []
+    result_pka_file_bonds = open(result_pka_pkafile[:-4] + ".bonds", "r")
+    for l in result_pka_file_bonds:
+        bonds.append(l.split())
+    result_pka_file_bonds.close()
+    return(bonds)
+
+
+def createdirs():
+    if platform.system() == 'Windows':
+        Newdir = os.getcwd() + "\Results_propka\\"
+    if platform.system() == 'Linux':
+        Newdir = os.getcwd() + "/Results_propka/"
+    if not os.path.exists(Newdir):
+        os.makedirs(Newdir)
+    return(Newdir)
+
+
+def openfiles(Newdir, filename, logtime, source):
+    if source == "upload":
+        result_pka_pkafile = filename.replace(".pdb", ".pka")
+        result_pka_input_pkafile = filename.replace(".pdb", ".propka_input")
+        result_log_name = "%s_Results.log" % (Newdir)
+        result_pka_file_stripped_name = filename.replace(".pdb", ".stripped")
+        result_pka_file_bonds_name = filename.replace(".pdb", ".bonds")
+    if source == "ID":
+        result_pka_pkafile = "%s%s%s.pka" % (Newdir, filename, logtime)
+        result_pka_input_pkafile = "%s%s%s.propka_input" % (Newdir, filename, logtime)
+        result_log_name = "%s_Results.log" % (Newdir)
+        result_pka_file_stripped_name = "%s%s%s.stripped" % (Newdir, filename, logtime)
+        result_pka_file_bonds_name = "%s%s%s.bonds" % (Newdir, filename, logtime)
+    if platform.system() == 'Windows':
+        filepath = "\\"
+    if platform.system() == 'Linux':
+        filepath = "/"
+    # Open the files
+    result_pka_file = open(result_pka_pkafile, "w")
+    result_input_pka_file = open(result_pka_input_pkafile, "w")
+    result_log = open(result_log_name, "a")
+    result_pka_file_stripped = open(result_pka_file_stripped_name, "w")
+    result_pka_file_bonds = open(result_pka_file_bonds_name, "w")
+    return(result_pka_file, result_input_pka_file, result_log, filepath, result_pka_pkafile, result_pka_file_stripped, result_pka_file_bonds)
+
+
+def ResiRange(resi):
+    resi = resi.split('.')
+    resiList = []
+    for i in resi:
+        if '-' in i:
+            tmp = i.split('-')
+            resiList.extend(list(range(int(tmp[0]), int(tmp[-1]) + 1)))
+        if '-' not in i:
+            resiList.append(int(i))
+    return(resiList)
+
+
+def ResnRange(resn):
+    resn_split = resn.split('.')
+    resn_range = [resnr.upper() for resnr in resn_split]
+    return(resn_range)
+
+
+def ChainRange(chain):
+    chainstring = chain.replace(".", "+").upper()
+    return(chainstring)
+
+
+def writepymolcmd(newmolecule, pkafile, verbose, makebonds):
+    list_results, ligands_results = importpropkaresult(pkafile)
+    # Now find the available bonding partners that pymol knows of
+    bonding_partners = []
+    bonding_partners_str = cmd.get_pdbstr("%s and resn * and not resn ASP+GLU+ARG+LYS+HIS+CYS+TYR+GLN+ASN+SER+THR+GLY+PHE+LEU+ALA+ILE+TRP+MET+PRO+VAL+HOH" % (newmolecule))
+    for i in range(len(bonding_partners_str.splitlines()) - 1):
+        bonding_partners_split = bonding_partners_str.splitlines()[i].split()
+        if bonding_partners_split[0] == "HETATM" or bonding_partners_split[0] == "ATOM":
+            bonding_partners_single = bonding_partners_split[3]
+            bonding_partners.append(bonding_partners_single)
+    bonding_partners = uniqifi(bonding_partners)
+    if verbose == 'yes':
+        print("And other possible bonding partners is: %s" % bonding_partners)
+    # Read in the bond file, if it exists
+    writebonds = "no"
+    if os.path.isfile(pkafile[:-4] + ".bonds") and makebonds == "yes":
+        bonds = importpropkabonds(pkafile)
+        writebonds = "yes"
+    # Open the pymol command file for writing
+    files_pka_pymol = openpymolfiles(pkafile)
+    result_pka_pymol = files_pka_pymol[0]
+    result_pka_pymol_name = files_pka_pymol[1]
+    # Make some dictionary for propka->pymol name conversion
+    dictio = {'ASP': 'CG', 'GLU': 'CD', 'ARG': 'CZ', 'LYS': 'NZ', 'HIS': 'CG', 'CYS': 'SG', 'TYR': 'OH', 'C-': 'C', 'N+': 'N', 'NTR': 'N', 'CTR': 'C', 'GLN': 'CD', 'ASN': 'CG', 'SER': 'OG', 'THR': 'OG1', 'GLY': 'CA', 'PHE': 'CZ', 'LEU': 'CG', 'ALA': 'CB', 'ILE': 'CD1', 'TRP': 'NE1', 'MET': 'SD', 'PRO': 'CG', 'VAL': 'CB'}
+    dictio2 = {'ASP': 'D', 'GLU': 'E', 'ARG': 'R', 'LYS': 'K', 'HIS': 'H', 'CYS': 'C', 'TYR': 'Y', 'C-': 'C-', 'N+': 'N+'}
+    # This list is from: http://en.wikipedia.org/wiki/Protein_pKa_calculations
+    pkaaminoacid = ['ASP', 'GLU', 'ARG', 'LYS', 'HIS', 'CYS', 'TYR']
+    pkadictio = {'ASP': 3.9, 'GLU': 4.3, 'ARG': 12.0, 'LYS': 10.5, 'HIS': 6.0, 'CYS': 8.3, 'TYR': 10.1}
+    # Now start write to the file.
+    # Try to make silent
+    # result_pka_pymol.write("cmd.feedback('disable','all','actions')\n")
+    # result_pka_pymol.write("cmd.feedback('disable','all','results')\n")
+    # Change the GUI width, to make the long names possible.
+    result_pka_pymol.write("cmd.set('internal_gui_width','360')\n")
+    # Set fonts
+    result_pka_pymol.write("cmd.set('label_font_id','12')\n")
+    result_pka_pymol.write("cmd.set('label_size','-0.5')\n")
+    result_pka_pymol.write("cmd.set('label_color','grey')\n")
+    # No auto zoom the new objects
+    result_pka_pymol.write("cmd.set('auto_zoom','off')\n")
+    # The name for the molecules are defined here
+    pkamolecule = "%spKa" % (newmolecule)
+    pkalabelmolecule = "%sLab" % (newmolecule)
+    # Create the groups now, so they come in order. They will be empty
+    result_pka_pymol.write("cmd.group('%sResi','Res*')\n" % (newmolecule))
+    result_pka_pymol.write("cmd.group('%sLigands','Lig*')\n" % (newmolecule))
+    if writebonds == "yes":
+        result_pka_pymol.write("cmd.group('%sBonds','%sBond*')\n" % (newmolecule, newmolecule))
+    # Create new empty pymol pka molecules. For pka atoms and its label. This is a "bucket" we where we will put in the atoms together.
+    result_pka_pymol.write("cmd.create('%s','None')\n" % (pkamolecule))
+    result_pka_pymol.write("cmd.create('%s','None')\n" % (pkalabelmolecule))
+    # Now make the pka atoms and alter, color and such
+    for l in list_results:
+        name = dictio[l[0]]
+        resn = dictio2[l[0]]
+        resi = l[1]
+        chain = l[2]
+        pka = l[3]
+        buried = l[4]
+        if "*" in pka:
+            pka = pka.replace("*", "")
+            comment = "*Coupled residue"
+        else:
+            comment = ""
+        if l[0] in pkaaminoacid:
+            pkadiff = (float(pka) - pkadictio[l[0]])
+            pkadiff = "(%s)" % pkadiff
+            if pka == "99.99":
+                pkadiff = ""
+        else:
+            pkadiff = ""
+        # Make the selection for which atom to copy
+        newselection = ("/%s//%s/%s/%s" % (newmolecule, chain, resi, name))
+        protselect = ("%sRes_%s%s%s" % (newmolecule, chain, resn, resi))
+        result_pka_pymol.write("cmd.select('%s','byres %s')\n" % (protselect, newselection))
+        result_pka_pymol.write("cmd.show('sticks','byres %s')\n" % (protselect))
+        # The temporary name
+        tempname = ("%s%s%s%s" % (pkamolecule, chain, resi, name))
+        tempnamelabel = ("%s%s%s%s" % (pkalabelmolecule, chain, resi, name))
+        tempselect = ("/%s//%s/%s" % (tempname, chain, resi))
+        tempselectlabel = ("/%s//%s/%s" % (tempnamelabel, chain, resi))
+        # Copy the atom, call it by the residue name
+        result_pka_pymol.write("cmd.create('%s','%s',quiet=1)\n" % (tempname, newselection))
+        # Alter the name and the b value of the newly created atom
+        result_pka_pymol.write("cmd.alter('%s','b=%s')\n" % (tempselect, pka))
+        result_pka_pymol.write("cmd.alter('%s','vdw=0.5')\n" % (tempselect))
+        result_pka_pymol.write("cmd.alter('%s','name=%s%s%s')\n" % (tempselect, '"', pka, '"'))
+        # Now create a fake label atom, and translate it
+        result_pka_pymol.write("cmd.create('%s','%s',quiet=1)\n" % (tempnamelabel, tempname))
+        movelabelxyz = (1.5, 0, 0)
+        result_pka_pymol.write("cmd.translate('[%s,%s,%s]','%s',camera=0)\n" % (movelabelxyz[0], movelabelxyz[1], movelabelxyz[2], tempnamelabel))
+        # Labelling alternate positions are not allowed, so we delete that attribute for the label atoms.
+        result_pka_pymol.write("cmd.alter('%s','alt=%s%s')\n" % (tempselectlabel, '"', '"'))
+        result_pka_pymol.write("cmd.label('%s','text_type=%spka=%s%s Bu:%s%s%s%s')\n" % (tempselectlabel, '"', pka, pkadiff, buried, '%', comment, '"'))
+        # Now put the atoms into a bucket of atoms
+        result_pka_pymol.write("cmd.create('%s','%s or (%s)',quiet=1)\n" % (pkamolecule, pkamolecule, tempselect))
+        result_pka_pymol.write("cmd.create('%s','%s or (%s)',quiet=1)\n" % (pkalabelmolecule, pkalabelmolecule, tempselectlabel))
+        # Remove the temporary atoms
+        result_pka_pymol.write("cmd.remove('%s')\n" % (tempname))
+        result_pka_pymol.write("cmd.remove('%s')\n" % (tempnamelabel))
+        # Delete the temporary molecule/selection
+        result_pka_pymol.write("cmd.delete('%s')\n" % (tempname))
+        result_pka_pymol.write("cmd.delete('%s')\n" % (tempnamelabel))
+    # Group the resi together
+    result_pka_pymol.write("cmd.group('%sResi','%sRes*')\n" % (newmolecule, newmolecule))
+    for l in ligands_results:
+        resn = l[0]
+        atom = l[1]
+        chain = l[2]
+        pka = l[3]
+        buried = l[4]
+        if verbose == 'yes':
+            print("Ligand. resn:%s atom:%s chain:%s pka:%s buried:%s" % (resn, atom, chain, pka, buried))
+        if Check_bonding_partners(bonding_partners, resn)[0]:
+            if "*" in pka:
+                pka = pka.replace("*", "")
+                comment = "*Coupled residue"
+            else:
+                comment = ""
+            # Make the selection for which atom to copy
+            ligselection = ("/%s and chain %s and resn %s and name %s" % (newmolecule, chain, resn, atom))
+            ligselect = ("%sLig_%s%s%s" % (newmolecule, chain, resn, atom))
+            result_pka_pymol.write("cmd.select('%s','%s')\n" % (ligselect, ligselection))
+            result_pka_pymol.write("cmd.show('sticks','byres %s')\n" % (ligselect))
+            result_pka_pymol.write("cmd.util.cbap('byres %s')\n" % (ligselect))
+            # The temporary name
+            tempname = ("%s%s%s%s" % (pkamolecule, chain, resn, atom))
+            tempnamelabel = ("%s%s%s%s" % (pkalabelmolecule, chain, resn, atom))
+            tempselect = ("/%s and chain %s and resn %s" % (tempname, chain, resn))
+            tempselectlabel = ("/%s and chain %s and resn %s" % (tempnamelabel, chain, resn))
+            # Copy the atom, call it by the residue name
+            result_pka_pymol.write("cmd.create('%s','%s',quiet=1)\n" % (tempname, ligselection))
+            # Alter the name and the b value of the newly created atom
+            result_pka_pymol.write("cmd.alter('%s','b=%s')\n" % (tempselect, pka))
+            result_pka_pymol.write("cmd.alter('%s','vdw=0.5')\n" % (tempselect))
+            result_pka_pymol.write("cmd.alter('%s','name=%s%s%s')\n" % (tempselect, '"', pka, '"'))
+            # Now create a fake label atom, and translate it
+            result_pka_pymol.write("cmd.create('%s','%s',quiet=1)\n" % (tempnamelabel, tempname))
+            movelabelxyz = (1.5, 0, 0)
+            result_pka_pymol.write("cmd.translate('[%s,%s,%s]','%s',camera=0)\n" % (movelabelxyz[0], movelabelxyz[1], movelabelxyz[2], tempnamelabel))
+            # Labelling alternate positions are not allowed, so we delete that attribute for the label atoms.
+            result_pka_pymol.write("cmd.alter('%s','alt=%s%s')\n" % (tempselectlabel, '"', '"'))
+            result_pka_pymol.write("cmd.label('%s','text_type=%spka=%s Bu:%s%s%s%s')\n" % (tempselectlabel, '"', pka, buried, '%', comment, '"'))
+            # Now put the atoms into a bucket of atoms
+            result_pka_pymol.write("cmd.create('%s','%s or (%s)',quiet=1)\n" % (pkamolecule, pkamolecule, tempselect))
+            result_pka_pymol.write("cmd.create('%s','%s or (%s)',quiet=1)\n" % (pkalabelmolecule, pkalabelmolecule, tempselectlabel))
+            # Remove the temporary atoms
+            result_pka_pymol.write("cmd.remove('%s')\n" % (tempname))
+            result_pka_pymol.write("cmd.remove('%s')\n" % (tempnamelabel))
+            # Delete the temporary molecule/selection
+            result_pka_pymol.write("cmd.delete('%s')\n" % (tempname))
+            result_pka_pymol.write("cmd.delete('%s')\n" % (tempnamelabel))
+    # Group the resi together
+    result_pka_pymol.write("cmd.group('%sLigands','%sLig*')\n" % (newmolecule, newmolecule))
+    # Finish the pka atoms, and show spheres
+    result_pka_pymol.write("cmd.show('spheres','%s')\n" % (pkamolecule))
+    result_pka_pymol.write("cmd.spectrum('b','red_white_blue',selection='%s',minimum='0',maximum='14')\n" % (pkamolecule))
+    result_pka_pymol.write("cmd.alter('%s and name 99.9','vdw=0.8')\n" % (pkamolecule))
+    result_pka_pymol.write("cmd.show('spheres','%s and name 99.9')\n" % (pkamolecule))
+    result_pka_pymol.write("cmd.color('sulfur','%s and name 99.9')\n" % (pkamolecule))
+    # Now we make the bonds
+    if writebonds == "yes":
+        Bondgroups = []
+        naturalaminoacids = ['ASP', 'GLU', 'ARG', 'LYS', 'HIS', 'CYS', 'TYR', 'NTR', 'N+', 'CTR', 'C-', 'GLN', 'ASN', 'SER', 'THR', 'GLY', 'PHE', 'LEU', 'ALA', 'ILE', 'TRP', 'MET', 'PRO', 'VAL']
+        for l in bonds:
+            if l[0] in naturalaminoacids:
+                name = dictio[l[0]]
+                resi = l[1]
+                chain = l[2]
+                desolvation = l[6][12:]
+                pkachange = l[11]
+                NBresi = l[8][3:]
+                NBchain = l[9]
+                NBbond = l[-1][:2]
+                if l[8][:3] in naturalaminoacids:
+                    NBname, cutoff = BondTypeName(dictio[l[8][:3]], NBbond)
+                    fromselection = ("/%s//%s/%s/%s" % (newmolecule, chain, resi, name))
+                    toselection = ("/%s//%s/%s/%s" % (newmolecule, NBchain, NBresi, NBname))
+                    if l[8][:3] == 'NTR':
+                        extind = cmd.identify("chain %s and name N" % (NBchain))[0]
+                        toselection = ("/%s and id %s and name N" % (newmolecule, extind))
+                        NBresi = "N+"
+                    if l[8][:3] == 'CTR':
+                        extind = cmd.identify("chain %s and name C" % (NBchain))[-1]
+                        toselection = ("/%s and id %s and name C" % (newmolecule, extind))
+                        NBresi = "C-"
+                    distname = ("%s_%s%s%s%s_%s_%s" % (newmolecule, chain, resi, NBchain, NBresi, NBbond, pkachange))
+                    result_pka_pymol.write("cmd.distance('%s','%s','%s'%s)\n" % (distname, fromselection, toselection, cutoff))
+                    result_pka_pymol.write("cmd.color('%s','%s')\n" % (SetDashColor(NBbond), distname))
+                    # result_pka_pymol.write("cmd.disable('%s')\n"%(distname))
+                    Bondgroups.append("%s%s" % (chain, resi))
+                if l[8][:3] not in naturalaminoacids and Check_bonding_partners(bonding_partners, l[8])[0]:
+                    cutoff = ""
+                    NBresn = Check_bonding_partners(bonding_partners, l[8])[1]
+                    NBname = l[8][len(NBresn):] + "*"
+                    fromselection = ("/%s//%s/%s/%s" % (newmolecule, chain, resi, name))
+                    toselection = ("/%s and chain %s and resn %s and name %s" % (newmolecule, NBchain, NBresn, NBname))
+                    if verbose == 'yes':
+                        print("Res->Ligand: (%s) -> (%s)" % (fromselection, toselection))
+                    result_pka_pymol.write("cmd.show('sticks','%s')\n" % (toselection))
+                    distname = ("%s_%s%s%s_%s_%s" % (newmolecule, chain, resi, NBresn, NBbond, pkachange))
+                    result_pka_pymol.write("cmd.distance('%s','%s','%s'%s)\n" % (distname, fromselection, toselection, cutoff))
+                    result_pka_pymol.write("cmd.color('%s','%s')\n" % (SetDashColor(NBbond), distname))
+                    # result_pka_pymol.write("cmd.disable('%s')\n"%(distname))
+                    Bondgroups.append("%s%s" % (chain, resi))
+            if l[0] in bonding_partners:
+                resn = l[0]
+                atom = l[1]
+                chain = l[2]
+                desolvation = l[6][12:]
+                pkachange = l[11]
+                NBresi = l[8][3:]
+                NBchain = l[9]
+                NBbond = l[-1][:2]
+                if not Check_bonding_partners(bonding_partners, l[8])[0]:
+                    NBname, cutoff = BondTypeName(dictio[l[8][:3]], NBbond)
+                    fromselection = ("/%s and chain %s and resn %s and name %s" % (newmolecule, chain, resn, atom))
+                    toselection = ("/%s//%s/%s/%s" % (newmolecule, NBchain, NBresi, NBname))
+                    if l[8][:3] == 'NTR':
+                        extind = cmd.identify("chain %s and name N" % (NBchain))[0]
+                        toselection = ("/%s and id %s and name N" % (newmolecule, extind))
+                        NBresi = "N+"
+                    if l[8][:3] == 'CTR':
+                        extind = cmd.identify("chain %s and name C" % (NBchain))[-1]
+                        toselection = ("/%s and id %s and name C" % (newmolecule, extind))
+                        NBresi = "C-"
+                    distname = ("%s_%s%s%s%s%s_%s_%s" % (newmolecule, chain, resn, atom, NBchain, NBresi, NBbond, pkachange))
+                    result_pka_pymol.write("cmd.distance('%s','%s','%s'%s)\n" % (distname, fromselection, toselection, cutoff))
+                    result_pka_pymol.write("cmd.color('%s','%s')\n" % (SetDashColor(NBbond), distname))
+                    Bondgroups.append("%s%s%s" % (chain, resn, atom))
+                    # result_pka_pymol.write("cmd.disable('%s')\n"%(distname))
+                if Check_bonding_partners(bonding_partners, l[8])[0]:
+                    cutoff = ""
+                    NBresn = Check_bonding_partners(bonding_partners, l[8])[1]
+                    NBname = l[8][len(NBresn):] + "*"
+                    fromselection = ("/%s and chain %s and resn %s and name %s" % (newmolecule, chain, resn, atom))
+                    toselection = ("/%s and chain %s and resn %s and name %s" % (newmolecule, NBchain, NBresn, NBname))
+                    if verbose == 'yes':
+                        print("Ligand->Ligand: (%s) -> (%s)" % (fromselection, toselection))
+                    result_pka_pymol.write("cmd.show('sticks','%s')\n" % (toselection))
+                    distname = ("%s_%s%s%s%s_%s_%s" % (newmolecule, chain, resn, atom, NBresn, NBbond, pkachange))
+                    result_pka_pymol.write("cmd.distance('%s','%s','%s'%s)\n" % (distname, fromselection, toselection, cutoff))
+                    result_pka_pymol.write("cmd.color('%s','%s')\n" % (SetDashColor(NBbond), distname))
+                    # result_pka_pymol.write("cmd.disable('%s')\n"%(distname))
+                    Bondgroups.append("%s%s%s" % (chain, resn, atom))
+        Bondgroups = uniqifi(Bondgroups)
+        for l in Bondgroups:
+            result_pka_pymol.write("cmd.group('%sBonds_%s','%s_%s*')\n" % (newmolecule, l, newmolecule, l))
+            result_pka_pymol.write("cmd.disable('%sBonds_%s')\n" % (newmolecule, l))
+        result_pka_pymol.write("cmd.group('%sBonds','%sBonds_*')\n" % (newmolecule, newmolecule))
+    result_pka_pymol.write("cmd.set('auto_zoom','on')\n")
+    # result_pka_pymol.write("cmd.feedback('enable','all','actions')\n")
+    # result_pka_pymol.write("cmd.feedback('enable','all','results')\n")
+    result_pka_pymol.close()
+    return(result_pka_pymol_name)
+
+
+def replace_all(text, dic):
+    for i, j in dic.items():
+        text = text.replace(i, j)
+    return(text)
+
+
+def uniqifi(seq, idfun=None):
+    # Order preserving
+    if idfun is None:
+        def idfun(x):
+            return x
+    seen = {}
+    result = []
+    for item in seq:
+        marker = idfun(item)
+        if marker in seen:
+            continue
+        seen[marker] = 1
+        result.append(item)
+    return(result)
+
+
+def BondTypeName(NBname, NBbond):
+    if NBbond == "SH":
+        cutoff = ""
+        return(NBname, cutoff)
+    if NBbond == "BH":
+        cutoff = ""
+        return("N", cutoff)
+    else:
+        cutoff = ""
+        return(NBname, cutoff)
+
+
+def Check_bonding_partners(bonding_partners, NBname):
+    answer = False
+    for l in bonding_partners:
+        if l in NBname:
+            answer = True
+            NBname = l
+            break
+        else:
+            answer = False
+    return(answer, NBname)
+
+
+def SetDashColor(NBbond):
+    if NBbond == "SH":
+        color = "brightorange"
+    if NBbond == "BH":
+        color = "lightorange"
+    if NBbond == "CC":
+        color = "red"
+    return(color)

.pymol/startup/show_contact_map.py

@@ -0,0 +1,58 @@
+#!/usr/bin/env python
+
+import tkFileDialog
+import tkSimpleDialog
+from pymol import cmd
+import math
+from distutils.util import strtobool
+
+def __init__(self):
+	self.menuBar.addmenuitem('Plugin', 'command', 'Contact map', label='Show contact map', command=lambda s=self: load_map(s))
+
+def load_map(app):
+	contacts_file = tkFileDialog.askopenfile(initialdir=".", title="Load contact map")
+	if contacts_file is None:
+		return
+	show_contact_map(file=contacts_file)
+
+def show_contact_map(selection="all", file=None, as_bonds=1, quiet=0):
+	""" Visualize contact map
+	"""
+	as_bonds = bool(as_bonds)
+	quiet = bool(quiet)
+	contacts = read_contacts(file)
+
+	for contact in contacts:	
+		name = "contact{}-{}".format(contact[0], contact[1])
+		seleA = "id %s and %s" % (contact[0], selection)
+		seleB = "id %s and %s" % (contact[1], selection)
+		
+		# skip non-exsistant
+		if cmd.select("seleA", seleA) == "0" or cmd.select("seleB", seleB) == "0":
+			print("%s to %s not found!" % (seleA, seleB))
+
+		if not quiet:
+			print("New bond: residue %s to %s" % (contact[0], contact[1]))
+
+		if as_bonds:
+			cmd.bond(seleA, seleB)
+		else:
+			cmd.distance(name, seleA, seleB)
+
+	cmd.delete('seleA')
+	cmd.delete('seleB')
+
+
+def read_contacts(contacts_file, separator=' '):
+	contacts = []
+	if not isinstance(contacts_file, file):
+		contacts_file = open(contacts_file, 'r')
+	for line in contacts_file.readlines():
+		split = line.split(",")
+		resA = int(split[0])
+		resB = int(split[1])
+		contacts.append([resA, resB])
+	contacts_file.close()
+	return contacts
+	
+cmd.extend("contact_map", show_contact_map)

.pymol/startup/spectrum_states.py

@@ -0,0 +1,92 @@
+'''
+http://pymolwiki.org/index.php/spectrum_states
+
+(c) 2011 Takanori Nakane and Thomas Holder
+
+License: BSD-2-Clause
+'''
+
+from __future__ import print_function
+
+from pymol import cmd, CmdException
+
+
+def spectrum_states(selection='all', representations='cartoon ribbon',
+                    color_list='blue cyan green yellow orange red',
+                    first=1, last=0, quiet=1):
+    '''
+DESCRIPTION
+
+    Color each state in a multi-state object different.
+
+USAGE
+
+    spectrum_states [ selection [, representations [, color_list [, first [, last ]]]]]
+
+ARGUMENTS
+
+    selection = string: object names (works with complete objects only)
+    {default: all}
+
+    representations = string: space separated list of representations
+    {default: cartoon ribbon}
+
+    color_list = string: space separated list of colors {default: blue cyan
+    green yellow orange red}
+
+SEE ALSO
+
+    spectrum, spectrumany
+    '''
+    from math import floor, ceil
+
+    first, last, quiet = int(first), int(last), int(quiet)
+    colors = color_list.split()
+    if len(colors) < 2:
+        print(' Error: please provide at least 2 colors')
+        raise CmdException
+
+    colvec = [cmd.get_color_tuple(i) for i in colors]
+
+    # filter for valid <repr>_color settings
+    settings = []
+    for r in representations.split():
+        if r[-1] == 's':
+            r = r[:-1]
+        s = r + '_color'
+        if s in cmd.setting.name_list:
+            settings.append(s)
+        elif not quiet:
+            print(' Warning: no such setting:', s)
+
+    # object names only
+    selection = ' '.join(cmd.get_object_list('(' + selection + ')'))
+    if cmd.count_atoms(selection) == 0:
+        print(' Error: empty selection')
+        raise CmdException
+
+    if last < 1:
+        last = cmd.count_states(selection)
+
+    val_range = int(last - first + 1)
+    if val_range < 2:
+        print(' Error: no spectrum possible, need more than 1 state')
+        raise CmdException
+
+    for i in range(val_range):
+        p = float(i) / (val_range - 1) * (len(colvec) - 1)
+        p0, p1 = int(floor(p)), int(ceil(p))
+        ii = (p - p0)
+        col_list = [colvec[p1][j] * ii + colvec[p0][j] * (1.0 - ii) for j in range(3)]
+        col_name = '0x%02x%02x%02x' % (col_list[0] * 255, col_list[1] * 255, col_list[2] * 255)
+        for s in settings:
+            cmd.set(s, col_name, selection, state=i + first)
+
+cmd.extend('spectrum_states', spectrum_states)
+
+# tab-completion of arguments
+cmd.auto_arg[0]['spectrum_states'] = cmd.auto_arg[0]['disable']
+cmd.auto_arg[1]['spectrum_states'] = [cmd.auto_arg[0]['show'][0], 'representation', ' ']
+cmd.auto_arg[2]['spectrum_states'] = [cmd.auto_arg[0]['color'][0], 'color', ' ']
+
+# vi:expandtab:smarttab

.pymol/startup/tmalign.py

@@ -0,0 +1,257 @@
+'''
+(c) 2011-2012 Thomas Holder, MPI for Developmental Biology
+'''
+
+from __future__ import print_function
+
+__author__ = 'Thomas Holder'
+__version__ = '1.1'
+__license__ = 'BSD-2-Clause'
+
+from pymol import cmd, CmdException
+
+
+def save_pdb_without_ter(filename, selection, **kwargs):
+    '''
+DESCRIPTION
+
+Save PDB file without TER records. External applications like TMalign and
+DynDom stop reading PDB files at TER records, which might be undesired in
+case of missing loops.
+    '''
+    v = cmd.get_setting_boolean('pdb_use_ter_records')
+    if v:
+        cmd.unset('pdb_use_ter_records')
+    cmd.save(filename, selection, **kwargs)
+    if v:
+        cmd.set('pdb_use_ter_records')
+
+
+def alignwithanymethod(mobile, target, methods='align super cealign tmalign',
+                       async=1, quiet=1):
+    '''
+DESCRIPTION
+
+Align copies of mobile to target with several alignment methods
+
+ARGUMENTS
+
+mobile = string: atom selection
+
+target = string: atom selection
+
+methods = string: space separated list of PyMOL commands which take
+arguments "mobile" and "target" (in any order) {default: align super
+cealign tmalign}
+    '''
+    import threading
+    import time
+    methods = methods.split()
+    async, quiet = int(async), int(quiet)
+    mobile_obj = cmd.get_object_list('first (' + mobile + ')')[0]
+
+    def myalign(method):
+        newmobile = cmd.get_unused_name(mobile_obj + '_' + method)
+        cmd.create(newmobile, mobile_obj)
+        start = time.time()
+        cmd.do('%s mobile=%s in %s, target=%s' % (method, newmobile, mobile, target))
+        if not quiet:
+            print('Finished: %s (%.2f sec)' % (method, time.time() - start))
+
+    for method in methods:
+        if async:
+            t = threading.Thread(target=myalign, args=(method,))
+            t.setDaemon(1)
+            t.start()
+        else:
+            myalign(method)
+
+
+def tmalign(mobile, target, args='', exe='TMalign', ter=0, transform=1, object=None, quiet=0):
+    '''
+DESCRIPTION
+
+TMalign wrapper
+
+Reference: Y. Zhang and J. Skolnick, Nucl. Acids Res. 2005 33, 2302-9
+http://zhanglab.ccmb.med.umich.edu/TM-align/
+
+USAGE
+
+tmalign mobile, target [, args [, exe ]]
+
+ARGUMENTS
+
+mobile, target = string: atom selections
+
+args = string: Extra arguments like -d0 5 -L 100
+
+exe = string: Path to TMalign executable {default: TMalign}
+
+ter = 0/1: If ter=0, then ignore chain breaks because TMalign will stop
+at first TER record {default: 0}
+
+SEE ALSO
+
+tmscore, mmalign
+    '''
+    import subprocess
+    import tempfile
+    import os
+    import re
+
+    ter, quiet = int(ter), int(quiet)
+
+    mobile_filename = tempfile.mktemp('.pdb', 'mobile')
+    target_filename = tempfile.mktemp('.pdb', 'target')
+    matrix_filename = tempfile.mktemp('.txt', 'matrix')
+    mobile_ca_sele = '(%s) and (not hetatm) and name CA and alt +A' % (mobile)
+    target_ca_sele = '(%s) and (not hetatm) and name CA and alt +A' % (target)
+
+    if ter:
+        save = cmd.save
+    else:
+        save = save_pdb_without_ter
+    save(mobile_filename, mobile_ca_sele)
+    save(target_filename, target_ca_sele)
+
+    exe = cmd.exp_path(exe)
+    args = [exe, mobile_filename, target_filename, '-m', matrix_filename] + args.split()
+
+    try:
+        process = subprocess.Popen(args, stdout=subprocess.PIPE)
+        lines = process.stdout.readlines()
+    except OSError:
+        print('Cannot execute "%s", please provide full path to TMscore or TMalign executable' % (exe))
+        raise CmdException
+    finally:
+        os.remove(mobile_filename)
+        os.remove(target_filename)
+
+    # TMalign >= 2012/04/17
+    if os.path.exists(matrix_filename):
+        lines += open(matrix_filename).readlines()
+        os.remove(matrix_filename)
+
+    r = None
+    re_score = re.compile(r'TM-score\s*=\s*(\d*\.\d*)')
+    rowcount = 0
+    matrix = []
+    line_it = iter(lines)
+    alignment = []
+    for line in line_it:
+        if 4 >= rowcount > 0:
+            if rowcount >= 2:
+                a = list(map(float, line.split()))
+                matrix.extend(a[2:5])
+                matrix.append(a[1])
+            rowcount += 1
+        elif line.lower().startswith(' -------- rotation matrix'):
+            rowcount = 1
+        elif line.startswith('(":" denotes'):
+            alignment = [line_it.next().rstrip() for i in range(3)]
+        else:
+            match = re_score.search(line)
+            if match is not None:
+                r = float(match.group(1))
+        if not quiet:
+            print(line.rstrip())
+
+    if not quiet:
+        for i in range(0, len(alignment[0]) - 1, 78):
+            for line in alignment:
+                print(line[i:i + 78])
+            print('')
+
+    assert len(matrix) == 3 * 4
+    matrix.extend([0, 0, 0, 1])
+
+    if int(transform):
+        cmd.transform_selection('byobject (%s)' % (mobile), matrix, homogenous=1)
+
+    # alignment object
+    if object is not None:
+        mobile_idx, target_idx = [], []
+        space = {'mobile_idx': mobile_idx, 'target_idx': target_idx}
+        cmd.iterate(mobile_ca_sele, 'mobile_idx.append("%s`%d" % (model, index))', space=space)
+        cmd.iterate(target_ca_sele, 'target_idx.append("%s`%d" % (model, index))', space=space)
+        for i, aa in enumerate(alignment[0]):
+            if aa == '-':
+                mobile_idx.insert(i, None)
+        for i, aa in enumerate(alignment[2]):
+            if aa == '-':
+                target_idx.insert(i, None)
+        if (len(mobile_idx) == len(target_idx) == len(alignment[2])):
+            cmd.rms_cur(
+                ' '.join(idx for (idx, m) in zip(mobile_idx, alignment[1]) if m in ':.'),
+                ' '.join(idx for (idx, m) in zip(target_idx, alignment[1]) if m in ':.'),
+                cycles=0, matchmaker=4, object=object)
+        else:
+            print('Could not load alignment object')
+
+    if not quiet and r is not None:
+        print('Found in output TM-score = %.4f' % (r))
+
+    return r
+
+
+def tmscore(mobile, target, args='', exe='TMscore', quiet=0, **kwargs):
+    '''
+DESCRIPTION
+
+TMscore wrapper
+
+Reference: Yang Zhang and Jeffrey Skolnick, Proteins 2004 57: 702-710
+http://zhanglab.ccmb.med.umich.edu/TM-score/
+
+ARGUMENTS
+
+mobile, target = string: atom selections
+
+args = string: Extra arguments like -d 5
+
+exe = string: Path to TMscore executable {default: TMscore}
+
+ter = 0/1: If ter=0, then ignore chain breaks because TMscore will stop
+at first TER record {default: 0}
+
+SEE ALSO
+
+tmalign, mmalign
+    '''
+    kwargs.pop('_self', None)
+    return tmalign(mobile, target, args, exe, quiet=quiet, **kwargs)
+
+
+def mmalign(mobile, target, args='', exe='MMalign', ter=0, transform=1, quiet=0):
+    '''
+DESCRIPTION
+
+MMalign wrapper
+
+Reference: S. Mukherjee and Y. Zhang, Nucleic Acids Research 2009; 37: e83
+http://zhanglab.ccmb.med.umich.edu/MM-align/
+
+SEE ALSO
+
+tmalign, tmscore
+    '''
+    return tmalign(mobile, target, args, exe, ter, transform, quiet=quiet)
+
+# pymol commands
+cmd.extend('alignwithanymethod', alignwithanymethod)
+cmd.extend('tmalign', tmalign)
+cmd.extend('tmscore', tmscore)
+cmd.extend('mmalign', tmalign)
+
+# autocompletion
+cmd.auto_arg[0].update({
+    'tmalign': cmd.auto_arg[0]['align'],
+    'tmscore': cmd.auto_arg[0]['align'],
+    'mmalign': cmd.auto_arg[0]['align'],
+})
+cmd.auto_arg[1].update({
+    'tmalign': cmd.auto_arg[1]['align'],
+    'tmscore': cmd.auto_arg[1]['align'],
+    'mmalign': cmd.auto_arg[1]['align'],
+})

.pymol/startup/visualize_dca_scores.py

@@ -0,0 +1,87 @@
+#!/usr/bin/env python
+
+import tkFileDialog
+import tkSimpleDialog
+from pymol import cmd
+import math
+
+DEFAULT_CUTOFF = 10
+DEFAULT_FLOOR = 0.6
+DEFAULT_MAX_CONTACTS = -1
+
+def __init__(self):
+    self.menuBar.addmenuitem('Plugin', 'command', 'MSA Scores', label='MSA Scores', command=lambda s=self: load_scores_dialog(s))
+
+def load_scores_dialog(app):
+    scores_file = tkFileDialog.askopenfile(initialdir=".", title="Load scores")
+    if scores_file is None:
+        return
+    cutoff = tkSimpleDialog.askinteger("Cutoff", "Backbone cutoff", initialvalue=DEFAULT_CUTOFF, minvalue=1)
+    if cutoff is None:
+        return
+    floor = tkSimpleDialog.askfloat("Minimal Score", "Minimal Score", initialvalue=DEFAULT_FLOOR)
+    if floor is None:
+        return
+    max_contacts = tkSimpleDialog.askfloat("Limit contacts", "Maximal contacts to show", initialvalue=DEFAULT_MAX_CONTACTS)
+    if max_contacts is None:
+        return
+    show_scores(scores=scores_file, cutoff=cutoff, floor=floor, max_contacts=max_contacts)
+
+def show_scores(selection='all', scores="scores.csv", floor=DEFAULT_FLOOR, cutoff=DEFAULT_CUTOFF, max_contacts=None):
+    """ Visualize score csv 
+        selection: Selection for visualization
+        floor: minimal score for visualization
+        cutoff: minimal distance between amino acids
+        scores: scores file for visualization
+    """
+    cutoff = int(cutoff)
+    floor = float(floor)
+    scores = read_scores(scores)
+    
+    max_contacts = int(max_contacts)
+    if max_contacts is None or max_contacts == -1:
+        max_contacts = len(scores)
+
+    # filter for cutoff
+    scores = list(filter(lambda x: x[1] - x[0] >= cutoff, scores))
+
+    # filter for score floor
+    scores = list(filter(lambda x: x[2] >= floor, scores))
+
+    # show distances
+    counter = 0
+    for score in scores:    
+        if counter > max_contacts:
+            print("max contacts reached! (%s contacts not shown)" % (len(scores)-int(max_contacts)))
+            break
+        print("new msa contact between %s and %s" % (score[0], score[1]))
+        name = "msa{}-{}_{:.3f}".format(score[0], score[1], score[2])
+        seleA = "resid %s and n. CA and %s" % (score[0], selection)
+        seleB = "resid %s and n. CA and %s" % (score[1], selection)
+        if cmd.select("seleA", seleA) == "0" or cmd.select("seleB", seleB) == "0":
+            continue
+        cmd.distance(name, seleA, seleB)
+        counter += 1
+        
+    cmd.delete('seleA')
+    cmd.delete('seleB')
+
+
+def read_scores(scores_file):
+    scores = []
+    if not isinstance(scores_file, file):
+        scores_file = open(scores_file, 'r')
+    for line in scores_file.readlines():
+        split = line.split(",")
+        resA = int(split[0])
+        resB = int(split[1])
+        score = float(split[2])
+        dist = abs(resB - resA)
+        scores.append([resA, resB, score, dist])
+    scores_file.close()
+    return scores
+
+def remove_scores():
+    cmd.delete("msa*")
+cmd.extend("remove_scores", remove_scores)    
+cmd.extend("show_scores", show_scores)

.pymolrc

@@ -0,0 +1,4 @@
+set opaque_background, 0
+log_open /tmp/pymol_log.pml
+set cartoon_side_chain_helper, 1
+set cartoon_highlight_color, grey50