Package rdkit :: Package Chem :: Package EState :: Module EState
[hide private]
[frames] | no frames]

Source Code for Module rdkit.Chem.EState.EState

  1  # $Id$ 
  2  # 
  3  # Copyright (C) 2002-2006 greg Landrum and Rational Discovery LLC 
  4  # 
  5  #   @@ All Rights Reserved @@ 
  6  #  This file is part of the RDKit. 
  7  #  The contents are covered by the terms of the BSD license 
  8  #  which is included in the file license.txt, found at the root 
  9  #  of the RDKit source tree. 
 10  # 
 11  """ Basic EState definitions 
 12   
 13  """ 
 14  from __future__ import print_function 
 15  import numpy 
 16  from rdkit import Chem 
 17   
 18   
19 -def GetPrincipleQuantumNumber(atNum):
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
37 -def EStateIndices(mol, force=True):
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
77 -def MaxEStateIndex(mol, force=1):
78 return max(EStateIndices(mol, force))
79 80 81 MaxEStateIndex.version = "1.0.0" 82 83
84 -def MinEStateIndex(mol, force=1):
85 return min(EStateIndices(mol, force))
86 87 88 MinEStateIndex.version = "1.0.0" 89 90
91 -def MaxAbsEStateIndex(mol, force=1):
92 return max([abs(x) for x in EStateIndices(mol, force)])
93 94 95 MaxAbsEStateIndex.version = "1.0.0" 96 97
98 -def MinAbsEStateIndex(mol, force=1):
99 return min([abs(x) for x in EStateIndices(mol, force)])
100 101 102 MinAbsEStateIndex.version = "1.0.0" 103 104
105 -def _exampleCode():
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__': # pragma: nocover 116 _exampleCode() 117