← Back to branch summary

~ubuntu-branches/debian/stretch/openbabel/stretch

« back to all changes in this revision

Viewing changes to src/chiral.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2008-07-22 23:54:58 UTC
  • mfrom: (3.1.10 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080722235458-3o606czluviz4akx
Tags: 2.2.0-2
* Upload to unstable.
* debian/control: Updated descriptions.
* debian/patches/gauss_cube_format.patch: New patch, makes the 
  gaussian cube format available again.
* debian/rules (DEB_DH_MAKESHLIBS_ARGS_libopenbabel3): Removed.
* debian/rules (DEB_CONFIGURE_EXTRA_FLAGS): Likewise.
* debian/libopenbabel3.install: Adjust formats directory to single 
  version hierarchy.

Show diffs side-by-side

added added

removed removed

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