Mercurial > repos > guerler > springsuite
comparison spring_package/Energy.py @ 39:172398348efd draft
"planemo upload commit 26b4018c88041ee0ca7c2976e0a012015173d7b6-dirty"
author | guerler |
---|---|
date | Fri, 22 Jan 2021 15:50:27 +0000 |
parents | 0be0af9e695d |
children |
comparison
equal
deleted
inserted
replaced
38:80a4b98121b6 | 39:172398348efd |
---|---|
1 import math | 1 import math |
2 from os.path import dirname, realpath | |
2 | 3 |
3 NTYPE = 21 | 4 NTYPE = 21 |
4 NDIST = 20 | 5 NDIST = 20 |
5 NSCALE = 2.0 | 6 NSCALE = 2.0 |
6 | 7 |
7 | 8 |
8 class Energy: | 9 class Energy: |
9 def __init__(self): | 10 def __init__(self): |
10 self.dfire = list() | 11 self.dfire = list() |
11 with open("spring_package/dfire/dfire.txt") as file: | 12 dirPath = dirname(realpath(__file__)) |
13 with open("%s/Energy.data" % dirPath) as file: | |
12 for line in file: | 14 for line in file: |
13 self.dfire.append(float(line)) | 15 self.dfire.append(float(line)) |
14 | 16 |
15 def get(self, moleculeA, moleculeB): | 17 def get(self, residuesA, residuesB): |
16 result = 0 | 18 result = 0 |
17 chainA = list(moleculeA.calpha.keys())[0] | 19 for atomA in residuesA: |
18 chainB = list(moleculeB.calpha.keys())[0] | 20 indexA = self.toResCode(atomA["alignedResidue"]) |
19 for i in moleculeA.calpha[chainA]: | 21 for atomB in residuesB: |
20 atomA = moleculeA.calpha[chainA][i] | 22 indexB = self.toResCode(atomB["alignedResidue"]) |
21 indexA = self.toResCode(atomA["residue"]) | 23 dist2 = ((atomA["x"] - atomB["x"]) ** 2 + |
22 for j in moleculeB.calpha[chainB]: | 24 (atomA["y"] - atomB["y"]) ** 2 + |
23 atomB = moleculeB.calpha[chainB][j] | 25 (atomA["z"] - atomB["z"]) ** 2) |
24 indexB = self.toResCode(atomB["residue"]) | |
25 dist2 = (atomA["x"] - atomB["x"]) ** 2 + (atomA["y"] - atomB["y"]) ** 2 + (atomA["z"] - atomB["z"]) ** 2 | |
26 dist = int((math.sqrt(dist2) * NSCALE)) | 26 dist = int((math.sqrt(dist2) * NSCALE)) |
27 if dist < NDIST: | 27 if dist < NDIST: |
28 index = indexA * NTYPE * NDIST + indexB * NDIST + dist | 28 index = indexA * NTYPE * NDIST + indexB * NDIST + dist |
29 result = result + self.dfire[index] | 29 result = result + self.dfire[index] |
30 return result | 30 return result |
31 | 31 |
32 def getClashes(self, moleculeA, moleculeB, minDist=5.0): | |
33 minDist = minDist ** 2 | |
34 clashes = 0 | |
35 chainA = list(moleculeA.calpha.keys())[0] | |
36 chainB = list(moleculeB.calpha.keys())[0] | |
37 calphaA = moleculeA.calpha[chainA] | |
38 calphaB = moleculeB.calpha[chainB] | |
39 lenA = len(calphaA.keys()) | |
40 lenB = len(calphaB.keys()) | |
41 if lenA > lenB: | |
42 temp = calphaB | |
43 calphaB = calphaA | |
44 calphaA = temp | |
45 lenA = len(calphaA.keys()) | |
46 lenB = len(calphaB.keys()) | |
47 for i in calphaA: | |
48 atomA = calphaA[i] | |
49 for j in calphaB: | |
50 atomB = calphaB[j] | |
51 dist2 = ((atomA["x"] - atomB["x"]) ** 2 + | |
52 (atomA["y"] - atomB["y"]) ** 2 + | |
53 (atomA["z"] - atomB["z"]) ** 2) | |
54 if dist2 < minDist: | |
55 clashes = clashes + 1 | |
56 break | |
57 return clashes / float(lenA) | |
58 | |
32 def toResCode(self, seq): | 59 def toResCode(self, seq): |
33 code = dict(ALA=0, CYS=1, ASP=2, GLU=3, PHE=4, GLY=5, HIS=6, ILE=7, LYS=8, LEU=9, MET=10, | 60 code = dict(A=0, C=1, D=2, E=3, F=4, G=5, H=6, I=7, K=8, L=9, M=10, |
34 ASN=11, PRO=12, GLN=13, ARG=14, SER=15, THR=16, VAL=17, TRP=18, TYR=19) | 61 N=11, P=12, Q=13, R=14, S=15, T=16, V=17, W=18, Y=19) |
35 return code[seq] if seq in code else 20 | 62 return code[seq] if seq in code else 20 |