← Back to branch summary

~vcs-imports/openbabel/trunk

« back to all changes in this revision

Viewing changes to chiral.cpp

  • Committer: ghutchis
  • Date: 2001-11-27 18:50:36 UTC
  • Revision ID: svn-v4:86f38270-e075-4da6-bf68-7dcd545c500b:openbabel/trunk:46
Initial revision

Show diffs side-by-side

added added

removed removed

Lines of Context:
chiral.cpp
 
1
/**********************************************************************
 
2
Copyright (C) 1998-2001 by OpenEye Scientific Software, Inc.
 
3
 
 
4
This program is free software; you can redistribute it and/or modify
 
5
it under the terms of the GNU General Public License as published by
 
6
the Free Software Foundation version 2 of the License.
 
7
 
 
8
This program is distributed in the hope that it will be useful,
 
9
but WITHOUT ANY WARRANTY; without even the implied warranty of
 
10
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
11
GNU General Public License for more details.
 
12
***********************************************************************/
 
13
 
 
14
#include <list>
 
15
#include "mol.h"
 
16
#include "oeutil.h"
 
17
 
 
18
using namespace OpenEye;
 
19
 
 
20
#include "matrix.h"
 
21
#include "chiral.h"
 
22
 
 
23
void OEMol::FindChiralCenters()
 
24
{
 
25
  if (HasChiralityPerceived()) return;
 
26
  SetChiralityPerceived();
 
27
 
 
28
  //do quick test to see if there are any possible chiral centers
 
29
  bool mayHaveChiralCenter=false;
 
30
  OEAtom *atom,*nbr;
 
31
  vector<OEAtom*>::iterator i;
 
32
  for (atom = BeginAtom(i);atom;atom = NextAtom(i))
 
33
    if (atom->GetHyb() == 3 && atom->GetHvyValence() >= 3)
 
34
      {
 
35
        mayHaveChiralCenter=true;
 
36
        break;
 
37
      }
 
38
 
 
39
  if (!mayHaveChiralCenter) return;
 
40
 
 
41
  OEBond *bond;
 
42
  vector<OEBond*>::iterator j;
 
43
  for (bond = BeginBond(j);bond;bond = NextBond(j))
 
44
    if (bond->IsWedge() || bond->IsHash())
 
45
        (bond->GetBeginAtom())->SetChiral();
 
46
 
 
47
#define INTMETHOD
 
48
 
 
49
#ifdef INTMETHOD
 
50
        vector<unsigned int> vgid;
 
51
        GetGIDVector(vgid);
 
52
        vector<unsigned int> tlist;
 
53
        vector<unsigned int>::iterator k;
 
54
#else //use Golender floating point method
 
55
        vector<float> gp;  
 
56
        GraphPotentials(*this,gp);
 
57
        vector<float> tlist;
 
58
        vector<float>::iterator k;
 
59
#endif
 
60
 
 
61
  bool ischiral;
 
62
  for (atom = BeginAtom(i);atom;atom = NextAtom(i))
 
63
    if (atom->GetHyb() == 3 && atom->GetHvyValence() >= 3 && !atom->IsChiral())
 
64
      {
 
65
                tlist.clear();
 
66
                ischiral = true;
 
67
 
 
68
                for (nbr = atom->BeginNbrAtom(j);nbr;nbr = atom->NextNbrAtom(j))
 
69
                {
 
70
                        for (k = tlist.begin();k != tlist.end();k++)
 
71
#ifdef INTMETHOD
 
72
                                if (vgid[nbr->GetIdx()-1] == *k)
 
73
#else
 
74
                                if (fabs(gp[nbr->GetIdx()]-*k) < 0.001)
 
75
#endif
 
76
                                        ischiral = false;
 
77
 
 
78
#ifdef INTMETHOD
 
79
                        if (ischiral) tlist.push_back(vgid[nbr->GetIdx()-1]);
 
80
#else
 
81
                        if (ischiral) tlist.push_back(gp[nbr->GetIdx()]);
 
82
#endif
 
83
                        else break;
 
84
                }
 
85
 
 
86
        if (ischiral) atom->SetChiral();
 
87
      }
 
88
}
 
89
 
 
90
void GetChirality(OEMol &mol, vector<int> &chirality)
 
91
{
 
92
  chirality.resize(mol.NumAtoms()+1);
 
93
  fill(chirality.begin(),chirality.end(),0);
 
94
 
 
95
  OEAtom *atom;
 
96
  vector<OEAtom*>::iterator i;
 
97
  for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i))
 
98
    if (atom->IsChiral())
 
99
      {
 
100
        float sv = CalcSignedVolume(mol,atom);
 
101
        if (sv < 0.0f)  chirality[atom->GetIdx()-1] = -1;
 
102
        else if (sv > 0.0)  chirality[atom->GetIdx()-1] = 1;
 
103
    }
 
104
}
 
105
 
 
106
// Calculate the signed volume for an atom.  If the atom has a valence of 3
 
107
// the coordinates of an attached hydrogen are calculated
 
108
 
 
109
float CalcSignedVolume(OEMol &mol,OEAtom *atm)
 
110
{
 
111
  Vector tmp_crd;
 
112
  vector<int> nbr_atms;
 
113
  vector<Vector> nbr_crds;
 
114
  float hbrad = etab.CorrectedBondRad(1,0);
 
115
        
 
116
  if (atm->GetHvyValence() < 3)
 
117
    {
 
118
      cerr << "Cannot calculate a signed volume for an atom with a heavy atom valence of " << atm->GetHvyValence() << endl;
 
119
      exit(0);
 
120
    }
 
121
 
 
122
  // Create a vector with the coordinates of the neighbor atoms
 
123
  OEAtom *nbr;
 
124
  vector<OEBond*>::iterator bint;
 
125
  for (nbr = atm->BeginNbrAtom(bint);nbr;nbr = atm->NextNbrAtom(bint))
 
126
    {
 
127
      nbr_atms.push_back(nbr->GetIdx());                
 
128
    }
 
129
  // sort the neighbor atoms to insure a consistent ordering
 
130
  sort(nbr_atms.begin(),nbr_atms.end());
 
131
  for (unsigned int i = 0; i < nbr_atms.size(); i++)
 
132
    {
 
133
      OEAtom *tmp_atm = mol.GetAtom(nbr_atms[i]);
 
134
      nbr_crds.push_back(tmp_atm->GetVector());
 
135
    }   
 
136
 
 
137
  // If we have three heavy atoms we need to calculate the position of the fourth
 
138
  if (atm->GetHvyValence() == 3)
 
139
    {           
 
140
      float bondlen = hbrad+etab.CorrectedBondRad(atm->GetAtomicNum(),atm->GetHyb());
 
141
      atm->GetNewBondVector(tmp_crd,bondlen);
 
142
      nbr_crds.push_back(tmp_crd);
 
143
    }
 
144
 
 
145
  return(signed_volume(nbr_crds[0],nbr_crds[1],nbr_crds[2],nbr_crds[3]));
 
146
}
 
147
 
 
148
// calculate a signed volume given a set of 4 coordinates
 
149
float signed_volume(const Vector &a, const Vector &b, const Vector &c, const Vector &d)
 
150
{
 
151
  Vector A,B,C;
 
152
  A = b-a;
 
153
  B = c-a;
 
154
  C = d-a;
 
155
  Matrix3x3 m(A,B,C);
 
156
  return(m.determinant());
 
157
}
 
158
 
 
159
// Calculate the Graph Potentials of a molecule
 
160
// based on 
 
161
// V.E. and Rozenblit, A.B. Golender
 
162
// Logical and Combinatorial Algorithms for Drug Design
 
163
// for an example see
 
164
// Walters, W. P., Yalkowsky, S. H., JCICS, 1996, 36(5), 1015-1017
 
165
 
 
166
void GraphPotentials(OEMol &mol, vector<float> &pot)
 
167
{
 
168
  float det;
 
169
 
 
170
  vector<vector<float> > g,c,h;
 
171
  construct_g_matrix(mol,g);
 
172
  invert_matrix(g,det); 
 
173
  construct_c_matrix(mol,c);
 
174
  mult_matrix(h,g,c);
 
175
  pot.resize(mol.NumAtoms()+1);
 
176
 
 
177
  for (unsigned int i = 0; i < mol.NumAtoms();i++) pot[i+1] = h[i][0];
 
178
}
 
179
 
 
180
 
 
181
// Construct the matrix G, which puts each atoms valence+1
 
182
// on the diagonal and and -1 on the off diagonal if two
 
183
// atoms are connected.
 
184
 
 
185
void construct_g_matrix(OEMol &mol, vector<vector<float> > &m)
 
186
{
 
187
  unsigned int i,j;
 
188
 
 
189
  OEAtom *atm1,*atm2;
 
190
  vector<OEAtom*>::iterator aint,bint;
 
191
 
 
192
  m.resize(mol.NumAtoms());
 
193
  for (i = 0; i < m.size(); i++)
 
194
    m[i].resize(mol.NumAtoms());
 
195
 
 
196
  for (atm1 = mol.BeginAtom(aint),i=0;atm1;atm1 = mol.NextAtom(aint),i++)
 
197
    for (atm2 = mol.BeginAtom(bint),j=0;atm2;atm2 = mol.NextAtom(bint),j++)
 
198
      {
 
199
        if (i == j)
 
200
          {
 
201
            m[i][j] = atm1->GetValence() + 1;
 
202
            m[i][j] += (float)atm1->GetAtomicNum()/10.0;
 
203
            m[i][j] += (float)atm1->GetHyb()/100.0;
 
204
          }
 
205
        else
 
206
          {
 
207
            if (atm1->IsConnected(atm2))
 
208
              m[i][j] = -1;
 
209
            else
 
210
              m[i][j] = 0;
 
211
          }
 
212
      }
 
213
}
 
214
 
 
215
// Construct the matrix C, which is simply a column vector
 
216
// consisting of the valence for each atom
 
217
void construct_c_matrix(OEMol &mol,vector<vector<float > > &m)
 
218
{
 
219
  unsigned int i;
 
220
  OEAtom *atm1;
 
221
  vector<OEAtom*>::iterator aint;
 
222
  
 
223
  m.resize(mol.NumAtoms());
 
224
  for (i = 0; i < m.size(); i++)
 
225
    m[i].resize(1);
 
226
  for (atm1 = mol.BeginAtom(aint),i=0;atm1;atm1 = mol.NextAtom(aint),i++)
 
227
    {
 
228
      m[i][0] = atm1->GetValence();
 
229
    }
 
230
}
 
231
 
 
232
 
 
233
 
 
234