ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-3.0/src/openbabel/chiral.cpp
Revision: 2440
Committed: Wed Nov 16 19:42:11 2005 UTC (18 years, 9 months ago) by tim
File size: 13711 byte(s)
Log Message:
adding openbabel

File Contents

# User Rev Content
1 tim 2440 /*********************************************************************
2     chiral.cpp - Detect chiral atoms and molecules.
3    
4     Copyright (C) 1998-2001 by OpenEye Scientific Software, Inc.
5     Some portions Copyright (C) 2001-2005 by Geoffrey R. Hutchison
6    
7     This file is part of the Open Babel project.
8     For more information, see <http://openbabel.sourceforge.net/>
9    
10     This program is free software; you can redistribute it and/or modify
11     it under the terms of the GNU General Public License as published by
12     the Free Software Foundation version 2 of the License.
13    
14     This program is distributed in the hope that it will be useful,
15     but WITHOUT ANY WARRANTY; without even the implied warranty of
16     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17     GNU General Public License for more details.
18     ***********************************************************************/
19    
20     #include <list>
21    
22     #include "mol.hpp"
23     #include "obutil.hpp"
24     #include "matrix.hpp"
25     #include "chiral.hpp"
26     #include "matrix3x3.hpp"
27    
28     using namespace std;
29     namespace OpenBabel
30     {
31    
32     //! Sets atom->IsChiral() to true for chiral atoms
33     void OBMol::FindChiralCenters()
34     {
35     if (HasChiralityPerceived())
36     return;
37     SetChiralityPerceived();
38    
39     obErrorLog.ThrowError(__FUNCTION__,
40     "Ran OpenBabel::FindChiralCenters", obAuditMsg);
41    
42     //do quick test to see if there are any possible chiral centers
43     bool mayHaveChiralCenter=false;
44     OBAtom *atom,*nbr;
45     vector<OBNodeBase*>::iterator i;
46     for (atom = BeginAtom(i);atom;atom = NextAtom(i))
47     if (atom->GetHyb() == 3 && atom->GetHvyValence() >= 3)
48     {
49     mayHaveChiralCenter=true;
50     break;
51     }
52    
53     if (!mayHaveChiralCenter)
54     return;
55    
56     OBBond *bond;
57     vector<OBEdgeBase*>::iterator j;
58     for (bond = BeginBond(j);bond;bond = NextBond(j))
59     if (bond->IsWedge() || bond->IsHash())
60     (bond->GetBeginAtom())->SetChiral();
61    
62     #define INTMETHOD
63    
64     #ifdef INTMETHOD
65    
66     vector<unsigned int> vgid;
67     GetGIDVector(vgid);
68     vector<unsigned int> tlist;
69     vector<unsigned int>::iterator k;
70     #else //use Golender floating point method
71    
72     vector<double> gp;
73     GraphPotentials(*this,gp);
74     vector<double> tlist;
75     vector<double>::iterator k;
76     #endif
77    
78     bool ischiral;
79     for (atom = BeginAtom(i);atom;atom = NextAtom(i))
80     if (atom->GetHyb() == 3 && atom->GetHvyValence() >= 3 && !atom->IsChiral())
81     {
82     tlist.clear();
83     ischiral = true;
84    
85     for (nbr = atom->BeginNbrAtom(j);nbr;nbr = atom->NextNbrAtom(j))
86     {
87     for (k = tlist.begin();k != tlist.end();k++)
88     #ifdef INTMETHOD
89    
90     if (vgid[nbr->GetIdx()-1] == *k)
91     #else
92    
93     if (fabs(gp[nbr->GetIdx()]-*k) < 0.001)
94     #endif
95    
96     ischiral = false;
97    
98     #ifdef INTMETHOD
99    
100     if (ischiral)
101     tlist.push_back(vgid[nbr->GetIdx()-1]);
102     #else
103    
104     if (ischiral)
105     tlist.push_back(gp[nbr->GetIdx()]);
106     #endif
107    
108     else
109     break;
110     }
111    
112     if (ischiral)
113     atom->SetChiral();
114     }
115     }
116    
117     // Seems to make a vector chirality become filled with array of +/- 1 for chiral atoms.
118     void GetChirality(OBMol &mol, std::vector<int> &chirality)
119     {
120     chirality.resize(mol.NumAtoms()+1);
121     fill(chirality.begin(),chirality.end(),0);
122    
123     OBAtom *atom;
124     vector<OBNodeBase*>::iterator i;
125     for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i))
126     if (atom->IsChiral())
127     {
128     if (!atom->HasChiralVolume())
129     {
130     double sv = CalcSignedVolume(mol,atom);
131     if (sv < 0.0)
132     {
133     chirality[atom->GetIdx()-1] = -1;
134     atom->SetNegativeStereo();
135     }
136     else if (sv > 0.0)
137     {
138     chirality[atom->GetIdx()-1] = 1;
139     atom->SetPositiveStereo();
140     }
141     }
142     else // already calculated signed volume (e.g., imported from somewhere)
143     {
144     if (atom ->IsPositiveStereo())
145     chirality[atom->GetIdx()-1] = 1;
146     else
147     chirality[atom->GetIdx()-1] = -1;
148     }
149     }
150     }
151    
152     int GetParity4Ref(vector<unsigned int> pref)
153     {
154     if(pref.size()!=4)return(-1); // should be given a vector of size 4.
155     int parity=0;
156     for (int i=0;i<3;i++) // do the bubble sort this many times
157     {
158     for(int j=0;j<3;j++) // iterate across the array 4th element has no
159     { // right hand neighbour so no need to sort
160     if (pref[j+1] < pref[j]) // compare the two neighbors
161     {
162     unsigned int tmp = pref[j]; // swap a[j] and a[j+1]
163     pref[j] = pref[j+1];
164     pref[j+1] = tmp;
165     parity++; // parity odd will invert stereochem
166     }
167     }
168     } // End Bubble Sort
169     return(parity%2);
170     }
171    
172     bool CorrectChirality(OBMol &mol, OBAtom *atm, atomreftype i, atomreftype o)
173     {
174     if (!atm->HasChiralitySpecified()) // if no chirality defined can't do any more for 0D
175     return(false);
176    
177     int parityI=0,parityO=0;
178     OBChiralData* cd=(OBChiralData*)atm->GetData(OBGenericDataType::ChiralData);
179     if ((cd->GetAtom4Refs(input)).size()!=4)return(false); // must have 4 refs
180     parityI=GetParity4Ref(cd->GetAtom4Refs(i)); // Gets Atom4Refs used to define the chirality
181     parityO=GetParity4Ref(cd->GetAtom4Refs(o));//GetsOutput parity.
182     /* switch (CHTYPE)
183     {
184     case SMILES: // SMILES always uses 1234 atom refs
185     parityO=0; // if parityO==parityI then clockwise in input = clockwise in output
186     break;
187     case MOLV3000: // MOLV3000 uses 1234 unless an H then 123H
188     if (atm->GetHvyValence()==3)
189     {
190     OBAtom *nbr;
191     int Hid=1000;// max Atom ID +1 should be used here
192     vector<int> nbr_atms;
193     vector<OBEdgeBase*>::iterator i;
194     for (nbr = atm->BeginNbrAtom(i);nbr;nbr = atm->NextNbrAtom(i))
195     {
196     if (nbr->IsHydrogen()){Hid=nbr->GetIdx();continue;}
197     nbr_atms.push_back(nbr->GetIdx());
198     }
199     sort(nbr_atms.begin(),nbr_atms.end());
200     nbr_atms.push_back(Hid);
201     int tmp[4];
202     for(int i=0;i<4;i++){tmp[i]=nbr_atms[i];}
203     parityO=GetParity4Ref(tmp);
204     }
205     else if (atm->GetHvyValence()==4)
206     parityO=0;
207     break;
208     default:
209     parityO=0;
210     }*/
211     if (parityO==parityI)
212     {//cout << "Parity is the same"<<endl;
213     return(true);
214     }
215     else if(parityO!=parityI) // Need to invert the Chirality which has been set
216     { //cout << "Parity is Opposite"<<endl;
217     if (atm->IsClockwise())
218     {atm->UnsetStereo();atm->SetAntiClockwiseStereo();}
219     else if (atm->IsAntiClockwise())
220     {atm->UnsetStereo();atm->SetClockwiseStereo();}
221     else
222     return(false);
223     return(true);
224     }
225     return false;
226     }
227    
228     //! Calculate the signed volume for an atom. If the atom has a valence of 3
229     //! the coordinates of an attached hydrogen are calculated
230     //! Puts attached Hydrogen last at the moment, like mol V3000 format.
231     double CalcSignedVolume(OBMol &mol,OBAtom *atm)
232     {
233     vector3 tmp_crd;
234     vector<unsigned int> nbr_atms;
235     vector<vector3> nbr_crds;
236     bool use_central_atom = false,is2D=false;
237     double hbrad = etab.CorrectedBondRad(1,0);
238    
239     if (!mol.Has3D()) //give peudo Z coords if mol is 2D
240     {
241     vector3 v,vz(0.0,0.0,1.0);
242     is2D = true;
243     OBAtom *nbr;
244     OBBond *bond;
245     vector<OBEdgeBase*>::iterator i;
246     for (bond = atm->BeginBond(i);bond;bond = atm->NextBond(i))
247     {
248     nbr = bond->GetEndAtom();
249     if (nbr != atm)
250     {
251     v = nbr->GetVector();
252     if (bond->IsWedge())
253     v += vz;
254     else
255     if (bond->IsHash())
256     v -= vz;
257    
258     nbr->SetVector(v);
259     }
260     else
261     {
262     nbr = bond->GetBeginAtom();
263     v = nbr->GetVector();
264     if (bond->IsWedge())
265     v -= vz;
266     else
267     if (bond->IsHash())
268     v += vz;
269    
270     nbr->SetVector(v);
271     }
272     }
273     }
274    
275     if (atm->GetHvyValence() < 3)
276     {
277     #ifdef HAVE_SSTREAM
278     stringstream errorMsg;
279     #else
280     strstream errorMsg;
281     #endif
282     errorMsg << "Cannot calculate a signed volume for an atom with a heavy atom valence of " << atm->GetHvyValence() << endl;
283     obErrorLog.ThrowError(__FUNCTION__, errorMsg.str(), obInfo);
284     return(0.0);
285     }
286    
287     // Create a vector with the coordinates of the neighbor atoms
288     // Also make a vector with Atom IDs
289     OBAtom *nbr;
290     vector<OBEdgeBase*>::iterator bint;
291     for (nbr = atm->BeginNbrAtom(bint);nbr;nbr = atm->NextNbrAtom(bint))
292     {
293     nbr_atms.push_back(nbr->GetIdx());
294     }
295     // sort the neighbor atoms to insure a consistent ordering
296     sort(nbr_atms.begin(),nbr_atms.end());
297     for (unsigned int i = 0; i < nbr_atms.size(); i++)
298     {
299     OBAtom *tmp_atm = mol.GetAtom(nbr_atms[i]);
300     nbr_crds.push_back(tmp_atm->GetVector());
301     }
302     /*
303     // If we have three heavy atoms we need to calculate the position of the fourth
304     if (atm->GetHvyValence() == 3)
305     {
306     double bondlen = hbrad+etab.CorrectedBondRad(atm->GetAtomicNum(),atm->GetHyb());
307     atm->GetNewBondVector(tmp_crd,bondlen);
308     nbr_crds.push_back(tmp_crd);
309     }
310     */
311     for(int j=0;j < nbr_crds.size();j++) // Checks for a neighbour having 0 co-ords (added hydrogen etc)
312     {
313     if (nbr_crds[j]==0 && use_central_atom==false)use_central_atom=true;
314     else if (nbr_crds[j]==0)
315     {
316     obErrorLog.ThrowError(__FUNCTION__, "More than 2 neighbours have 0 co-ords when attempting 3D chiral calculation", obInfo);
317     }
318     }
319    
320     // If we have three heavy atoms we can use the chiral center atom itself for the fourth
321     // will always give same sign (for tetrahedron), magnitude will be smaller.
322     if(nbr_atms.size()==3 || use_central_atom==true)
323     {
324     nbr_crds.push_back(atm->GetVector());
325     nbr_atms.push_back(mol.NumAtoms()+1); // meed to add largest number on end to work
326     }
327     OBChiralData* cd=(OBChiralData*)atm->GetData(OBGenericDataType::ChiralData); //Set the output atom4refs to the ones used
328     if(cd==NULL)
329     {
330     cd = new OBChiralData;
331     atm->SetData(cd);
332     }
333     cd->SetAtom4Refs(nbr_atms,calcvolume);
334    
335     //re-zero psuedo-coords
336     if (is2D)
337     {
338     vector3 v;
339     OBAtom *atom;
340     vector<OBNodeBase*>::iterator k;
341     for (atom = mol.BeginAtom(k);atom;atom = mol.NextAtom(k))
342     {
343     v = atom->GetVector();
344     v.SetZ(0.0);
345     atom->SetVector(v);
346     }
347     }
348    
349     return(signed_volume(nbr_crds[0],nbr_crds[1],nbr_crds[2],nbr_crds[3]));
350     }
351    
352     //! Calculate a signed volume given a set of 4 coordinates
353     double signed_volume(const vector3 &a, const vector3 &b, const vector3 &c, const vector3 &d)
354     {
355     vector3 A,B,C;
356     A = b-a;
357     B = c-a;
358     C = d-a;
359     matrix3x3 m(A,B,C);
360     return(m.determinant());
361     }
362    
363     //! \brief Calculate the Graph Potentials of a molecule
364     //!
365     //! based on
366     //! V.E. and Rozenblit, A.B. Golender
367     //! <em>Logical and Combinatorial Algorithms for Drug Design</em>. \n
368     //! For an example see:
369     //! Walters, W. P., Yalkowsky, S. H., \em JCICS, 1996, 36(5), 1015-1017.
370     //! <a href="http://dx.doi.org/10.1021/ci950278o">DOI: 10.1021/ci950278o</a>
371     void GraphPotentials(OBMol &mol, std::vector<double> &pot)
372     {
373     double det;
374    
375     vector<vector<double> > g,c,h;
376     construct_g_matrix(mol,g);
377     invert_matrix(g,det);
378     construct_c_matrix(mol,c);
379     mult_matrix(h,g,c);
380     pot.resize(mol.NumAtoms()+1);
381    
382     for (unsigned int i = 0; i < mol.NumAtoms();i++)
383     pot[i+1] = h[i][0];
384     }
385    
386    
387     //! Construct the matrix G, which puts each atoms valence+1
388     //! on the diagonal and and -1 on the off diagonal if two
389     //! atoms are connected.
390     void construct_g_matrix(OBMol &mol, std::vector<std::vector<double> > &m)
391     {
392     unsigned int i,j;
393    
394     OBAtom *atm1,*atm2;
395     vector<OBNodeBase*>::iterator aint,bint;
396    
397     m.resize(mol.NumAtoms());
398     for (i = 0; i < m.size(); i++)
399     m[i].resize(mol.NumAtoms());
400    
401     for (atm1 = mol.BeginAtom(aint),i=0;atm1;atm1 = mol.NextAtom(aint),i++)
402     for (atm2 = mol.BeginAtom(bint),j=0;atm2;atm2 = mol.NextAtom(bint),j++)
403     {
404     if (i == j)
405     {
406     m[i][j] = atm1->GetValence() + 1;
407     m[i][j] += (double)atm1->GetAtomicNum()/10.0;
408     m[i][j] += (double)atm1->GetHyb()/100.0;
409     }
410     else
411     {
412     if (atm1->IsConnected(atm2))
413     m[i][j] = -1;
414     else
415     m[i][j] = 0;
416     }
417     }
418     }
419    
420     //! Construct the matrix C, which is simply a column vector
421     //! consisting of the valence for each atom
422     void construct_c_matrix(OBMol &mol,std::vector<std::vector<double > > &m)
423     {
424     unsigned int i;
425     OBAtom *atm1;
426     vector<OBNodeBase*>::iterator aint;
427    
428     m.resize(mol.NumAtoms());
429     for (i = 0; i < m.size(); i++)
430     m[i].resize(1);
431     for (atm1 = mol.BeginAtom(aint),i=0;atm1;atm1 = mol.NextAtom(aint),i++)
432     {
433     m[i][0] = atm1->GetValence();
434     }
435     }
436    
437     } // namespace OpenBabel
438    
439     //! \file chiral.cpp
440     //! \brief Detect chiral atoms and molecules.