-
Notifications
You must be signed in to change notification settings - Fork 1
/
coreRMSD.py
52 lines (38 loc) · 1.19 KB
/
coreRMSD.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
import openbabel
from rdkit import Chem
from rdkit.Chem import AllChem
import math
def getRMSD(distances): #get final rmsd value
numDistances=len(distances)
itrDistance=0
deviations=0
while itrDistance<numDistances:
deviations=deviations+distances[itrDistance]
itrDistance=itrDistance+1
return round(math.sqrt(deviations/numDistances),3)
def getDistance(pos1,pos2): #calculate the distance square of two atoms
return (pos1.x-pos2.x)**2+(pos1.y-pos2.y)**2+(pos1.z-pos2.z)**2
obConversion=openbabel.OBConversion()
obConversion.SetInFormat("mol2")
print sb
obmol=openbabel.OBMol()
notatend=obConversion.ReadFile(obmol,"csr_ladi-results.mol2")
while notatend:
#rdmol=openbabel.OBMolToRWMol(obmol)
obmol=openbabel.OBMol()
notatend=obConversion.Read(obmol)
m1=Chem.MolFromMol2File("test1.mol2")
m2=Chem.MolFromMol2File("test2.mol2")
core=Chem.MolFromMol2File("core.mol2")
core1=Chem.ReplaceSidechains(m1,core)
core2=Chem.ReplaceSidechains(m2,core)
atoms1=core1.GetAtoms()
atoms2=core2.GetAtoms()
i=0
myarray=[]
while i < len(atoms1):
pos1=core1.GetConformer().GetAtomPosition(i)
pos2=core2.GetConformer().GetAtomPosition(i)
myarray.append(getDistance(pos1,pos2))
i=i+1
print getRMSD(myarray)