Source Code for Module Chem.EState.EState

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