1
2
3
4
5
6
7
8
9
10
11 """ Basic EState definitions
12
13 """
14 from __future__ import print_function
15 import numpy
16 from rdkit import Chem
17
18
20 """ Get principal quantum number for atom number """
21 if atNum <= 2:
22 return 1
23 elif atNum <= 10:
24 return 2
25 elif atNum <= 18:
26 return 3
27 elif atNum <= 36:
28 return 4
29 elif atNum <= 54:
30 return 5
31 elif atNum <= 86:
32 return 6
33 else:
34 return 7
35
36
38 """ returns a tuple of EState indices for the molecule
39
40 Reference: Hall, Mohney and Kier. JCICS _31_ 76-81 (1991)
41
42 """
43 if not force and hasattr(mol, '_eStateIndices'):
44 return mol._eStateIndices
45
46 tbl = Chem.GetPeriodicTable()
47 nAtoms = mol.GetNumAtoms()
48 Is = numpy.zeros(nAtoms, numpy.float)
49 for i in range(nAtoms):
50 at = mol.GetAtomWithIdx(i)
51 atNum = at.GetAtomicNum()
52 d = at.GetDegree()
53 if d > 0:
54 h = at.GetTotalNumHs()
55 dv = tbl.GetNOuterElecs(atNum) - h
56 N = GetPrincipleQuantumNumber(atNum)
57 Is[i] = (4. / (N * N) * dv + 1) / d
58 dists = Chem.GetDistanceMatrix(mol, useBO=0, useAtomWts=0)
59 dists += 1
60 accum = numpy.zeros(nAtoms, numpy.float)
61 for i in range(nAtoms):
62 for j in range(i + 1, nAtoms):
63 p = dists[i, j]
64 if p < 1e6:
65 tmp = (Is[i] - Is[j]) / (p * p)
66 accum[i] += tmp
67 accum[j] -= tmp
68
69 res = accum + Is
70 mol._eStateIndices = res
71 return res
72
73
74 EStateIndices.version = '1.0.0'
75
76
79
80
81 MaxEStateIndex.version = "1.0.0"
82
83
86
87
88 MinEStateIndex.version = "1.0.0"
89
90
93
94
95 MaxAbsEStateIndex.version = "1.0.0"
96
97
100
101
102 MinAbsEStateIndex.version = "1.0.0"
103
104
106 """ Example code for calculating E-state indices """
107 smis = ['CCCC', 'CCCCC', 'CCCCCC', 'CC(N)C(=O)O', 'CC(N)C(=O)[O-].[Na+]']
108 for smi in smis:
109 m = Chem.MolFromSmiles(smi)
110 print(smi)
111 inds = EStateIndices(m)
112 print('\t', inds)
113
114
115 if __name__ == '__main__':
116 _exampleCode()
117