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)20std::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