diff spring_package/Energy.py @ 39:172398348efd draft

"planemo upload commit 26b4018c88041ee0ca7c2976e0a012015173d7b6-dirty"

--- a/spring_package/Energy.py	Wed Nov 25 17:38:24 2020 +0000
+++ b/spring_package/Energy.py	Fri Jan 22 15:50:27 2021 +0000
@@ -1,4 +1,5 @@
 import math
+from os.path import dirname, realpath
 
 NTYPE = 21
 NDIST = 20
@@ -8,28 +9,54 @@
 class Energy:
     def __init__(self):
         self.dfire = list()
-        with open("spring_package/dfire/dfire.txt") as file:
+        dirPath = dirname(realpath(__file__))
+        with open("%s/Energy.data" % dirPath) as file:
             for line in file:
                 self.dfire.append(float(line))
 
-    def get(self, moleculeA, moleculeB):
+    def get(self, residuesA, residuesB):
         result = 0
-        chainA = list(moleculeA.calpha.keys())[0]
-        chainB = list(moleculeB.calpha.keys())[0]
-        for i in moleculeA.calpha[chainA]:
-            atomA = moleculeA.calpha[chainA][i]
-            indexA = self.toResCode(atomA["residue"])
-            for j in moleculeB.calpha[chainB]:
-                atomB = moleculeB.calpha[chainB][j]
-                indexB = self.toResCode(atomB["residue"])
-                dist2 = (atomA["x"] - atomB["x"]) ** 2 + (atomA["y"] - atomB["y"]) ** 2 + (atomA["z"] - atomB["z"]) ** 2
+        for atomA in residuesA:
+            indexA = self.toResCode(atomA["alignedResidue"])
+            for atomB in residuesB:
+                indexB = self.toResCode(atomB["alignedResidue"])
+                dist2 = ((atomA["x"] - atomB["x"]) ** 2 +
+                         (atomA["y"] - atomB["y"]) ** 2 +
+                         (atomA["z"] - atomB["z"]) ** 2)
                 dist = int((math.sqrt(dist2) * NSCALE))
                 if dist < NDIST:
                     index = indexA * NTYPE * NDIST + indexB * NDIST + dist
                     result = result + self.dfire[index]
         return result
 
+    def getClashes(self, moleculeA, moleculeB, minDist=5.0):
+        minDist = minDist ** 2
+        clashes = 0
+        chainA = list(moleculeA.calpha.keys())[0]
+        chainB = list(moleculeB.calpha.keys())[0]
+        calphaA = moleculeA.calpha[chainA]
+        calphaB = moleculeB.calpha[chainB]
+        lenA = len(calphaA.keys())
+        lenB = len(calphaB.keys())
+        if lenA > lenB:
+            temp = calphaB
+            calphaB = calphaA
+            calphaA = temp
+            lenA = len(calphaA.keys())
+            lenB = len(calphaB.keys())
+        for i in calphaA:
+            atomA = calphaA[i]
+            for j in calphaB:
+                atomB = calphaB[j]
+                dist2 = ((atomA["x"] - atomB["x"]) ** 2 +
+                         (atomA["y"] - atomB["y"]) ** 2 +
+                         (atomA["z"] - atomB["z"]) ** 2)
+                if dist2 < minDist:
+                    clashes = clashes + 1
+                    break
+        return clashes / float(lenA)
+
     def toResCode(self, seq):
-        code = dict(ALA=0, CYS=1, ASP=2, GLU=3, PHE=4, GLY=5, HIS=6, ILE=7, LYS=8, LEU=9, MET=10,
-                    ASN=11, PRO=12, GLN=13, ARG=14, SER=15, THR=16, VAL=17, TRP=18, TYR=19)
+        code = dict(A=0, C=1, D=2, E=3, F=4, G=5, H=6, I=7, K=8, L=9, M=10,
+                    N=11, P=12, Q=13, R=14, S=15, T=16, V=17, W=18, Y=19)
         return code[seq] if seq in code else 20