Source Code for Module rdkit.Chem.EState.EState

 1  # $Id: EState.py 1528 2010-09-26 17:04:37Z glandrum $ 
 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  import numpy 
15  from rdkit import Chem 
16   
17 -def GetPrincipleQuantumNumber(atNum): 
18    if atNum<=2: return 1 
19    elif atNum <= 10: return 2 
20    elif atNum <= 18: return 3 
21    elif atNum <= 36: return 4 
22    elif atNum <= 54: return 5 
23    elif atNum <= 86: return 6 
24    else: return 7 
25     
26   
27 -def EStateIndices(mol,force=1): 
28    """ returns a tuple of EState indices for the molecule 
29   
30      Reference: Hall, Mohney and Kier. JCICS _31_ 76-81 (1991) 
31   
32    """ 
33    if not force and hasattr(mol,'_eStateIndices'): 
34      return mol._eStateIndices 
35       
36    tbl = Chem.GetPeriodicTable() 
37    nAtoms = mol.GetNumAtoms() 
38    Is = numpy.zeros(nAtoms,numpy.float) 
39    for i in range(nAtoms): 
40      at = mol.GetAtomWithIdx(i) 
41      atNum = at.GetAtomicNum() 
42      d = at.GetDegree() 
43      if d>0: 
44        h = at.GetTotalNumHs() 
45        dv = tbl.GetNOuterElecs(atNum)-h 
46        N = GetPrincipleQuantumNumber(atNum) 
47        Is[i] = (4./(N*N) * dv + 1)/d 
48    dists = Chem.GetDistanceMatrix(mol,useBO=0,useAtomWts=0) 
49    dists += 1 
50    accum = numpy.zeros(nAtoms,numpy.float) 
51    for i in range(nAtoms): 
52      for j in range(i+1,nAtoms): 
53        p = dists[i,j] 
54        if p < 1e6: 
55          tmp = (Is[i]-Is[j])/(p*p) 
56          accum[i] += tmp 
57          accum[j] -= tmp 
58   
59    res = accum+Is 
60    mol._eStateIndices=res 
61    return res 
62  EStateIndices.version='1.0.0' 
63   
64  if __name__ =='__main__': 
65    smis = ['CCCC','CCCCC','CCCCCC','CC(N)C(=O)O','CC(N)C(=O)[O-].[Na+]'] 
66    for smi in smis: 
67      m = Chem.MolFromSmiles(smi) 
68      print smi 
69      inds = EStateIndices(m) 
70      print '\t',inds 
71