1 // $Id$
2 //
3 //  Copyright (C) 2013 Greg Landrum
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 
12 #include <GraphMol/RDKitBase.h>
13 #include <GraphMol/Descriptors/MolDescriptors.h>
14 #include <GraphMol/Descriptors/Lipinski.h>
15 #include <vector>
16 #include <algorithm>
17 
18 namespace RDKit {
19 namespace Descriptors {
calcMQNs(const ROMol & mol,bool force)20 std::vector<unsigned int> calcMQNs(const ROMol& mol, bool force) {
21   RDUNUSED_PARAM(force);
22   // FIX: use force value to enable caching
23   std::vector<unsigned int> res(42, 0);
24 
25   // ---------------------------------------------------
26   // atom-centered things
27   // Note: We're not doing exactly the same thing
28   //       as the original paper on polarity counts
29   //       since we're using different donor and acceptor
30   //       definitions.
31   ROMol::VERTEX_ITER atBegin, atEnd;
32   boost::tie(atBegin, atEnd) = mol.getVertices();
33   while (atBegin != atEnd) {
34     const Atom* at = mol[*atBegin];
35     ++atBegin;
36     unsigned int nHs = at->getTotalNumHs();
37     unsigned int nRings = mol.getRingInfo()->numAtomRings(at->getIdx());
38     switch (at->getAtomicNum()) {
39       case 0:
40       case 1:
41         break;
42       case 6:
43         res[0]++;
44         break;
45       case 9:
46         res[1]++;
47         break;
48       case 17:
49         res[2]++;
50         break;
51       case 35:
52         res[3]++;
53         break;
54       case 53:
55         res[4]++;
56         break;
57       case 16:
58         res[5]++;
59         break;
60       case 15:
61         res[6]++;
62         break;
63       case 7:
64         if (!nRings) {
65           res[7]++;
66         } else {
67           res[8]++;
68         }
69         if (at->getDegree() != 4) {
70           res[19]++;  // number of acceptor sites
71           res[20]++;  // number of acceptor atoms
72         }
73         if (nHs) {
74           res[21] += nHs;  // number of donor sites
75           res[22]++;       // number of donor atoms
76         }
77         break;
78       case 8:
79         if (!nRings) {
80           res[9]++;
81         } else {
82           res[10]++;
83         }
84         res[20]++;  // number of acceptor atoms
85         if (at->getFormalCharge() != -1) {
86           res[19] += 2;  // number of acceptor sites
87         } else {
88           res[19] += 3;  // number of acceptor sites
89         }
90         if (nHs) {
91           res[21] += nHs;  // number of donor sites
92           res[22]++;       // number of donor atoms
93         }
94         break;
95       default:
96         break;
97     }
98 
99     if (at->getFormalCharge() > 0) {
100       res[24]++;  // positive charges
101     } else if (at->getFormalCharge() < 0) {
102       res[23]++;  // negative charges
103     }
104 
105     if (at->getAtomicNum() != 1) {
106       switch (at->getDegree()) {
107         case 1:
108           res[25]++;
109           break;
110         case 2:
111           if (!nRings) {
112             res[26]++;
113           } else {
114             res[29]++;
115           }
116           break;
117         case 3:
118           if (!nRings) {
119             res[27]++;
120           } else {
121             res[30]++;
122           }
123           break;
124         case 4:
125           if (!nRings) {
126             res[28]++;
127           } else {
128             res[31]++;
129           }
130           break;
131       }
132       if (nRings >= 2) {
133         res[40]++;
134       }
135     }
136   }
137   res[11] = mol.getNumHeavyAtoms();
138 
139   // ---------------------------------------------------
140   // bond counts:
141   unsigned int nAromatic = 0;
142   ROMol::EDGE_ITER firstB, lastB;
143   boost::tie(firstB, lastB) = mol.getEdges();
144   while (firstB != lastB) {
145     const Bond* bond = mol[*firstB];
146     if (bond->getIsAromatic()) {
147       ++nAromatic;
148     }
149     unsigned int nRings = mol.getRingInfo()->numBondRings(bond->getIdx());
150     switch (bond->getBondType()) {
151       case Bond::SINGLE:
152         if (!nRings) {
153           res[12]++;
154         } else {
155           res[15]++;
156         }
157         break;
158       case Bond::DOUBLE:
159         if (!nRings) {
160           res[13]++;
161         } else {
162           res[16]++;
163         }
164         break;
165       case Bond::TRIPLE:
166         if (!nRings) {
167           res[14]++;
168         } else {
169           res[17]++;
170         }
171         break;
172       default:
173         break;
174     }
175     if (nRings >= 2) {
176       res[41]++;
177     }
178     ++firstB;
179   }
180   // rather than do the work to kekulize the molecule, we cheat
181   // by just dividing the number of aromatic bonds evenly among the
182   // cyclic single bond and cyclic double bond bins and give any
183   // remainder to the single bonds
184   res[15] += nAromatic / 2;
185   res[16] += nAromatic / 2;
186   if (nAromatic % 2) {
187     res[15]++;
188   }
189   res[18] = calcNumRotatableBonds(mol);
190 
191   // ---------------------------------------------------
192   //  ring size counts
193   for (const auto& iv : mol.getRingInfo()->atomRings()) {
194     if (iv.size() < 10) {
195       res[iv.size() + 29]++;
196     } else {
197       res[39]++;
198     }
199   }
200 
201   return res;
202 }
203 }  // end of namespace Descriptors
204 }  // end of namespace RDKit
205