import Numeric import LinearAlgebra import math from pymol import cmd #Must have a single atom as input def findAtomFromSelection(sel): model = cmd.get_model(sel) a = Numeric.array(model.atom[0].coord) return a def getTransformationMatrix(a_sel, b_sel, c_sel): a_orig = findAtomFromSelection(a_sel) b_orig = findAtomFromSelection(b_sel) c_orig = findAtomFromSelection(c_sel) a = Numeric.zeros(3, 'f') b = Numeric.zeros(3, 'f') c = Numeric.zeros(3, 'f') matrix = Numeric.zeros((4, 4), 'f') transformToOrigo(a_orig, b_orig, c_orig, a, b, c, matrix) matrix_t = LinearAlgebra.inverse(matrix) return matrix_t def listOfTransformationMatrix(matrix): list = [] list.append(matrix[0,0]) list.append(matrix[0,1]) list.append(matrix[0,2]) list.append(matrix[0,3]) list.append(matrix[1,0]) list.append(matrix[1,1]) list.append(matrix[1,2]) list.append(matrix[1,3]) list.append(matrix[2,0]) list.append(matrix[2,1]) list.append(matrix[2,2]) list.append(matrix[2,3]) list.append(matrix[3,0]) list.append(matrix[3,1]) list.append(matrix[3,2]) list.append(matrix[3,3]) return list #Computes euclidian distance between two points in three dimensions def euclidianDistance(a, b): try: assert len(a) is 3 assert len(b) is 3 return (math.sqrt((a[0]-b[0])**2 + (a[1]-b[1])**2 + (a[2]-b[2])**2)) except AssertionError: print 'euclidianDistance: a or b is not of three dimensions.' sys.exit(1) def cross(v, w): try: assert len(v) is 3 assert len(w) is 3 cross = Numeric.array([0.0, 0.0, 0.0]) cross[0] = v[1]*w[2] - v[2]*w[1] cross[1] = v[2]*w[0] - v[0]*w[2] cross[2] = v[0]*w[1] - v[1]*w[0] return cross except AssertionError: print 'cross: input vectors are not of dimension 3' sys.exit(1) #a, b and c-orig are coordinates, a, b and c coordinates, trans_matrix 4*4 float matrix #all of them Numeric.array def transformToOrigo(a_orig, b_orig, c_orig, a, b, c, trans_matrix): bx = 0.0 by = 0.0 #Vectors ab = Numeric.array([(b_orig[0]-a_orig[0]), (b_orig[1]-a_orig[1]), (b_orig[2]-a_orig[2])]) ac = Numeric.array([(c_orig[0]-a_orig[0]), (c_orig[1]-a_orig[1]), (c_orig[2]-a_orig[2])]) bc = Numeric.array([(c_orig[0]-b_orig[0]), (c_orig[1]-b_orig[1]), (c_orig[2]-b_orig[2])]) #Vector lengths ab_length = euclidianDistance(a_orig, b_orig) ac_length = euclidianDistance(a_orig, c_orig) bc_length = euclidianDistance(b_orig, c_orig) #Unit vectors x_unit = Numeric.zeros(3, 'f') y_unit = Numeric.zeros(3, 'f') z_unit = Numeric.zeros(3, 'f') #Final new coordinate for c c[0] = ac_length #Unit vector for x plane x_unit = ac/ac_length alpha = math.degrees(math.acos(Numeric.dot(ab, ac)/(ab_length*ac_length))) bx = ab_length*math.cos(math.radians(alpha)) by = math.sin(math.radians(alpha))*ab_length #Final new coordinates for b b[0] = bx b[1] = by #Unit vector for y plane y_unit = (ab - (bx*x_unit))/by #Unit vector for z plane z_unit = cross(x_unit, y_unit) trans_matrix[0][0] = x_unit[0] trans_matrix[1][0] = x_unit[1] trans_matrix[2][0] = x_unit[2] trans_matrix[3][0] = 0.0 trans_matrix[0][1] = y_unit[0] trans_matrix[1][1] = y_unit[1] trans_matrix[2][1] = y_unit[2] trans_matrix[3][1] = 0.0 trans_matrix[0][2] = z_unit[0] trans_matrix[1][2] = z_unit[1] trans_matrix[2][2] = z_unit[2] trans_matrix[3][2] = 0.0 trans_matrix[0][3] = a_orig[0] trans_matrix[1][3] = a_orig[1] trans_matrix[2][3] = a_orig[2] trans_matrix[3][3] = 1.0 def transformCoordinate(coord, trans_matrix): c = Numeric.array([coord[0], coord[1], coord[2], 1], 'f') return Numeric.matrixmultiply(trans_matrix, c)