Package rdkit :: Package Chem :: Package Features :: Module FeatDirUtilsRD

[hide private]
[frames] | no frames]

# Source Code for Module rdkit.Chem.Features.FeatDirUtilsRD

      1  # $Id$
      2  #
      3  # Copyright (C) 2004-2008 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  from rdkit import Geometry
     12  from rdkit import Chem
     13  import numpy
     14  import math
     15
     16  # BIG NOTE: we are going assume atom IDs starting from 0 instead of 1
     17  # for all the functions in this file. This is so that they
     18  # are reasonably indepedent of the combicode. However when using
     19  # with combicode the caller needs to make sure the atom IDs from combicode
     20  # are corrected before feeding them in here.
     21
     22 -def cross(v1,v2):
     23    res = numpy.array([ v1[1]*v2[2] - v1[2]*v2[1],
     24                  -v1[0]*v2[2] + v1[2]*v2[0],
     25                  v1[0]*v2[1] - v1[1]*v2[0]],numpy.double)
     26    return res
     27
     28 -def findNeighbors(atomId, adjMat):
     29    """
     30    Find the IDs of the neighboring atom IDs
     31
     32    ARGUMENTS:
     33    atomId - atom of interest
     34    adjMat - adjacency matrix for the compound
     35    """
     36    nbrs = []
     37    for i,nid in enumerate(adjMat[atomId]):
     38      if nid >= 1 :
     39        nbrs.append(i)
     40    return nbrs
     41
     42 -def _findAvgVec(conf, center, nbrs) :
     43    # find the average of the normalized vectors going from the center atoms to the
     44    # neighbors
     45    # the average vector is also normalized
     46    avgVec = 0
     47    for nbr in nbrs:
     48      nid = nbr.GetIdx()
     49      pt = conf.GetAtomPosition(nid)
     50      pt -= center
     51      pt.Normalize()
     52      if (avgVec == 0) :
     53        avgVec = pt
     54      else :
     55        avgVec += pt
     56
     57    avgVec.Normalize()
     58    return avgVec
     59
     60 -def GetAromaticFeatVects(conf, featAtoms, featLoc, scale=1.5):
     61    """
     62    Compute the direction vector for an aromatic feature
     63
     64    ARGUMENTS:
     65       conf - a conformer
     66       featAtoms - list of atom IDs that make up the feature
     67       featLoc - location of the aromatic feature specified as point3d
     68       scale - the size of the direction vector
     69    """
     70    dirType = 'linear'
     71    head = featLoc
     72    ats = [conf.GetAtomPosition(x) for x in featAtoms]
     73
     74    p0 = ats[0]
     75    p1 = ats[1]
     76    v1 = p0-head
     77    v2 = p1-head
     78    norm1 = v1.CrossProduct(v2)
     79    norm1.Normalize()
     80    norm1 *= scale
     81    #norm2 = norm1
     82    norm2 = head-norm1
     83    norm1 += head
     84    return ( (head,norm1),(head,norm2) ), dirType
     85
     86 -def ArbAxisRotation(theta,ax,pt):
     87    theta = math.pi*theta/180
     88    c = math.cos(theta)
     89    s = math.sin(theta)
     90    t = 1-c
     91    X = ax.x
     92    Y = ax.y
     93    Z = ax.z
     94    mat = [ [t*X*X+c, t*X*Y+s*Z, t*X*Z-s*Y],
     95            [t*X*Y-s*Z,t*Y*Y+c,t*Y*Z+s*X],
     96            [t*X*Z+s*Y,t*Y*Z-s*X,t*Z*Z+c] ]
     97    mat = numpy.array(mat)
     98    if isinstance(pt,Geometry.Point3D):
     99      pt = numpy.array((pt.x,pt.y,pt.z))
    100      tmp = numpy.dot(mat,pt)
    101      res=Geometry.Point3D(tmp[0],tmp[1],tmp[2])
    102    elif isinstance(pt,list) or isinstance(pt,tuple):
    103      pts = pt
    104      res = []
    105      for pt in pts:
    106        pt = numpy.array((pt.x,pt.y,pt.z))
    107        tmp = numpy.dot(mat,pt)
    108        res.append(Geometry.Point3D(tmp[0],tmp[1],tmp[2]))
    109    else:
    110      res=None
    111    return res
    112
    113
    114
    115
    116
    117 -def GetAcceptor2FeatVects(conf, featAtoms, scale=1.5):
    118    """
    119    Get the direction vectors for Acceptor of type 2
    120
    121    This is the acceptor with two adjacent heavy atoms. We will special case a few things here.
    122    If the acceptor atom is an oxygen we will assume a sp3 hybridization
    123    the acceptor directions (two of them)
    124    reflect that configurations. Otherwise the direction vector in plane with the neighboring
    125    heavy atoms
    126
    127    ARGUMENTS:
    128        featAtoms - list of atoms that are part of the feature
    129        scale - length of the direction vector
    130    """
    131    assert len(featAtoms) == 1
    132    aid = featAtoms[0]
    133    cpt = conf.GetAtomPosition(aid)
    134
    135    mol = conf.GetOwningMol()
    136    nbrs = list(mol.GetAtomWithIdx(aid).GetNeighbors())
    137    hydrogens = []
    138    tmp=[]
    139    while len(nbrs):
    140      nbr = nbrs.pop()
    141      if nbr.GetAtomicNum()==1:
    142        hydrogens.append(nbr)
    143      else:
    144        tmp.append(nbr)
    145    nbrs = tmp
    146    assert len(nbrs) == 2
    147
    148    bvec = _findAvgVec(conf, cpt, nbrs)
    149    bvec *= (-1.0*scale)
    150
    151    if (mol.GetAtomWithIdx(aid).GetAtomicNum()==8):
    152      # assume sp3
    153      # we will create two vectors by rotating bvec by half the tetrahedral angle in either directions
    154      v1 = conf.GetAtomPosition(nbrs[0].GetIdx())
    155      v1 -= cpt
    156      v2 = conf.GetAtomPosition(nbrs[1].GetIdx())
    157      v2 -= cpt
    158      rotAxis = v1 - v2
    159      rotAxis.Normalize()
    160      bv1 = ArbAxisRotation(54.5, rotAxis, bvec)
    161      bv1 += cpt
    162      bv2 = ArbAxisRotation(-54.5, rotAxis, bvec)
    163      bv2 += cpt
    164      return ((cpt, bv1), (cpt, bv2),), 'linear'
    165    else :
    166      bvec += cpt
    167      return ((cpt, bvec),), 'linear'
    168
    169 -def _GetTetrahedralFeatVect(conf,aid,scale):
    170    mol = conf.GetOwningMol()
    171
    172    cpt = conf.GetAtomPosition(aid)
    173    nbrs = mol.GetAtomWithIdx(aid).GetNeighbors()
    174    if not _checkPlanarity(conf,cpt,nbrs,tol=0.1):
    175      bvec = _findAvgVec(conf, cpt, nbrs)
    176      bvec *= (-1.0*scale)
    177      bvec += cpt
    178      res = ((cpt,bvec),)
    179    else:
    180      res = ()
    181    return res
    182
    183 -def GetDonor3FeatVects(conf, featAtoms, scale=1.5):
    184    """
    185    Get the direction vectors for Donor of type 3
    186
    187    This is a donor with three heavy atoms as neighbors. We will assume
    188    a tetrahedral arrangement of these neighbors. So the direction we are seeking
    189    is the last fourth arm of the sp3 arrangment
    190
    191    ARGUMENTS:
    192      featAtoms - list of atoms that are part of the feature
    193      scale - length of the direction vector
    194    """
    195    assert len(featAtoms) == 1
    196    aid = featAtoms[0]
    197
    198    tfv = _GetTetrahedralFeatVect(conf,aid,scale)
    199    return tfv, 'linear'
    200
    201 -def GetAcceptor3FeatVects(conf, featAtoms, scale=1.5):
    202    """
    203    Get the direction vectors for Donor of type 3
    204
    205    This is a donor with three heavy atoms as neighbors. We will assume
    206    a tetrahedral arrangement of these neighbors. So the direction we are seeking
    207    is the last fourth arm of the sp3 arrangment
    208
    209    ARGUMENTS:
    210      featAtoms - list of atoms that are part of the feature
    211      scale - length of the direction vector
    212    """
    213    assert len(featAtoms) == 1
    214    aid = featAtoms[0]
    215    tfv = _GetTetrahedralFeatVect(conf,aid,scale)
    216    return tfv, 'linear'
    217
    218 -def _findHydAtoms(nbrs, atomNames):
    219    hAtoms = []
    220    for nid in nbrs:
    221      if atomNames[nid] == 'H':
    222        hAtoms.append(nid)
    223    return hAtoms
    224
    225 -def _checkPlanarity(conf, cpt, nbrs, tol=1.0e-3):
    226    assert len(nbrs) == 3
    227    v1 = conf.GetAtomPosition(nbrs[0].GetIdx())
    228    v1 -= cpt
    229    v2 = conf.GetAtomPosition(nbrs[1].GetIdx())
    230    v2 -= cpt
    231    v3 = conf.GetAtomPosition(nbrs[2].GetIdx())
    232    v3 -= cpt
    233    normal = v1.CrossProduct(v2)
    234    dotP = abs(v3.DotProduct(normal))
    235    if (dotP <= tol) :
    236      return 1
    237    else :
    238      return 0
    239
    240 -def GetDonor2FeatVects(conf, featAtoms, scale=1.5) :
    241    """
    242    Get the direction vectors for Donor of type 2
    243
    244    This is a donor with two heavy atoms as neighbors. The atom may are may not have
    245    hydrogen on it. Here are the situations with the neighbors that will be considered here
    246      1. two heavy atoms and two hydrogens: we will assume a sp3 arrangement here
    247      2. two heavy atoms and one hydrogen: this can either be sp2 or sp3
    248      3. two heavy atoms and no hydrogens
    249
    250    ARGUMENTS:
    251      featAtoms - list of atoms that are part of the feature
    252      scale - length of the direction vector
    253    """
    254    assert len(featAtoms) == 1
    255    aid = featAtoms[0]
    256    mol = conf.GetOwningMol()
    257
    258    cpt = conf.GetAtomPosition(aid)
    259
    260    # find the two atoms that are neighbors of this atoms
    261    nbrs = list(mol.GetAtomWithIdx(aid).GetNeighbors())
    262    assert len(nbrs) >= 2
    263
    264    hydrogens = []
    265    tmp=[]
    266    while len(nbrs):
    267      nbr = nbrs.pop()
    268      if nbr.GetAtomicNum()==1:
    269        hydrogens.append(nbr)
    270      else:
    271        tmp.append(nbr)
    272    nbrs = tmp
    273
    274    if len(nbrs) == 2:
    275      # there should be no hydrogens in this case
    276      assert len(hydrogens) == 0
    277      # in this case the direction is the opposite of the average vector of the two neighbors
    278      bvec = _findAvgVec(conf, cpt, nbrs)
    279      bvec *= (-1.0*scale)
    280      bvec += cpt
    281      return ((cpt, bvec),), 'linear'
    282    elif len(nbrs) == 3:
    283      assert len(hydrogens) == 1
    284      # this is a little more tricky we have to check if the hydrogen is in the plane of the
    285      # two heavy atoms (i.e. sp2 arrangement) or out of plane (sp3 arrangement)
    286
    287      # one of the directions will be from hydrogen atom to the heavy atom
    288      hid = hydrogens[0].GetIdx()
    289      bvec = conf.GetAtomPosition(hid)
    290      bvec -= cpt
    291      bvec.Normalize()
    292      bvec *= scale
    293      bvec += cpt
    294      if _checkPlanarity(conf, cpt, nbrs):
    295        # only the hydrogen atom direction needs to be used
    296        return ((cpt, bvec),), 'linear'
    297      else :
    298        # we have a non-planar configuration - we will assume sp3 and compute a second direction vector
    299        ovec = _findAvgVec(conf, cpt, nbrs)
    300        ovec *= (-1.0*scale)
    301        ovec += cpt
    302        return ((cpt, bvec), (cpt, ovec),), 'linear'
    303
    304    elif len(nbrs) >= 4 :
    305      # in this case we should have two or more hydrogens we will simple use there directions
    306      res = []
    307      for hid in hydrogens:
    308        bvec = conf.GetAtomPosition(hid)
    309        bvec -= cpt
    310        bvec.Normalize()
    311        bvec *= scale
    312        bvec += cpt
    313        res.append((cpt, bvec))
    314      return tuple(res), 'linear'
    315
    316 -def GetDonor1FeatVects(conf, featAtoms, scale=1.5) :
    317    """
    318    Get the direction vectors for Donor of type 1
    319
    320    This is a donor with one heavy atom. It is not clear where we should we should be putting the
    321    direction vector for this. It should probably be a cone. In this case we will just use the
    322    direction vector from the donor atom to the heavy atom
    323
    324    ARGUMENTS:
    325
    326      featAtoms - list of atoms that are part of the feature
    327      scale - length of the direction vector
    328    """
    329    assert len(featAtoms) == 1
    330    aid = featAtoms[0]
    331    mol = conf.GetOwningMol()
    332    nbrs = mol.GetAtomWithIdx(aid).GetNeighbors()
    333
    334    # find the neighboring heavy atom
    335    hnbr = -1
    336    for nbr in nbrs:
    337      if nbr.GetAtomicNum()!=1:
    338        hnbr = nbr.GetIdx()
    339        break
    340
    341    cpt = conf.GetAtomPosition(aid)
    342    v1 = conf.GetAtomPosition(hnbr)
    343    v1 -= cpt
    344    v1.Normalize()
    345    v1 *= (-1.0*scale)
    346    v1 += cpt
    347    return ((cpt, v1),), 'cone'
    348
    349 -def GetAcceptor1FeatVects(conf, featAtoms, scale=1.5) :
    350    """
    351    Get the direction vectors for Acceptor of type 1
    352
    353    This is a acceptor with one heavy atom neighbor. There are two possibilities we will
    354    consider here
    355    1. The bond to the heavy atom is a single bond e.g. CO
    356       In this case we don't know the exact direction and we just use the inversion of this bond direction
    357       and mark this direction as a 'cone'
    358    2. The bond to the heavy atom is a double bond e.g. C=O
    359       In this case the we have two possible direction except in some special cases e.g. SO2
    360       where again we will use bond direction
    361
    362    ARGUMENTS:
    363      featAtoms - list of atoms that are part of the feature
    364      scale - length of the direction vector
    365    """
    366    assert len(featAtoms) == 1
    367    aid = featAtoms[0]
    368    mol = conf.GetOwningMol()
    369    nbrs = mol.GetAtomWithIdx(aid).GetNeighbors()
    370
    371    cpt = conf.GetAtomPosition(aid)
    372
    373    # find the adjacent heavy atom
    374    heavyAt = -1
    375    for nbr in nbrs:
    376      if nbr.GetAtomicNum()!=1:
    377        heavyAt = nbr
    378        break
    379
    380    singleBnd = mol.GetBondBetweenAtoms(aid,heavyAt.GetIdx()).GetBondType() > Chem.BondType.SINGLE
    381
    382    # special scale - if the heavy atom is a sulfur (we should proabably check phosphorous as well)
    383    sulfur = heavyAt.GetAtomicNum()==16
    384
    385    if singleBnd or sulfur:
    386      v1 = conf.GetAtomPosition(heavyAt.GetIdx())
    387      v1 -= cpt
    388      v1.Normalize()
    389      v1 *= (-1.0*scale)
    390      v1 += cpt
    391      return ((cpt, v1),), 'cone'
    392    else :
    393      # ok in this case we will assume that
    394      # heavy atom is sp2 hybridized and the direction vectors (two of them)
    395      # are in the same plane, we will find this plane by looking for one
    396      # of the neighbors of the heavy atom
    397      hvNbrs = heavyAt.GetNeighbors()
    398      hvNbr = -1
    399      for nbr in hvNbrs:
    400        if nbr.GetIdx() != aid:
    401          hvNbr = nbr
    402          break
    403
    404      pt1 = conf.GetAtomPosition(hvNbr.GetIdx())
    405      v1 = conf.GetAtomPosition(heavyAt.GetIdx())
    406      pt1 -= v1
    407      v1 -= cpt
    408      rotAxis = v1.CrossProduct(pt1)
    409      rotAxis.Normalize()
    410      bv1 = ArbAxisRotation(120, rotAxis, v1)
    411      bv1.Normalize()
    412      bv1 *= scale
    413      bv1 += cpt
    414      bv2 = ArbAxisRotation(-120, rotAxis, v1)
    415      bv2.Normalize()
    416      bv2 *= scale
    417      bv2 += cpt
    418      return ((cpt, bv1), (cpt, bv2),), 'linear'
    419