1 // $Id$
2 //
3 // Copyright (C) 2003-2011 Greg Landrum and Rational Discovery LLC
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 #include "GasteigerCharges.h"
12 #include <RDGeneral/types.h>
13 #include <GraphMol/RDKitBase.h>
14 #include <GraphMol/ROMol.h>
15 #include <GraphMol/MolOps.h>
16 #include "GasteigerParams.h"
17
18 namespace Gasteiger {
19 using namespace RDKit;
20 /*! \brief split the formal charge across atoms of same type if we have a
21 *conjugated system
22 *
23 * This function is called before the charge equivalization iteration is
24 *started for the
25 * Gasteiger charges. If any of the atom involved in conjugated system have
26 *formal charges
27 * set on them, this charges is equally distributed across atoms of the same
28 *type in that
29 * conjugated system. So for example the two nitrogens in the benzamidine
30 *system start the iteration
31 * with equal charges of 0.5
32 */
splitChargeConjugated(const ROMol & mol,DOUBLE_VECT & charges)33 void splitChargeConjugated(const ROMol &mol, DOUBLE_VECT &charges) {
34 int aix;
35 int natms = mol.getNumAtoms();
36 INT_VECT marker;
37 INT_VECT_CI mci;
38 int aax, yax;
39 double formal;
40 const Atom *at, *aat, *yat;
41 for (aix = 0; aix < natms; aix++) {
42 at = mol.getAtomWithIdx(aix);
43 formal = at->getFormalCharge();
44 // std::cout << aix << " formal charges:" << formal << "\n";
45 marker.resize(0);
46 if ((fabs(formal) > EPS_DOUBLE) && (fabs(charges[aix]) < EPS_DOUBLE)) {
47 marker.push_back(aix);
48 ROMol::OEDGE_ITER bnd1, end1, bnd2, end2;
49 boost::tie(bnd1, end1) = mol.getAtomBonds(at);
50 while (bnd1 != end1) {
51 if (mol[*bnd1]->getIsConjugated()) {
52 aax = mol[*bnd1]->getOtherAtomIdx(aix);
53 aat = mol.getAtomWithIdx(aax);
54 boost::tie(bnd2, end2) = mol.getAtomBonds(aat);
55 while (bnd2 != end2) {
56 if ((*bnd1) != (*bnd2)) {
57 if (mol[*bnd2]->getIsConjugated()) {
58 yax = mol[*bnd2]->getOtherAtomIdx(aax);
59 yat = mol.getAtomWithIdx(yax);
60 if (at->getAtomicNum() == yat->getAtomicNum()) {
61 formal += yat->getFormalCharge();
62 marker.push_back(yax);
63 }
64 }
65 }
66 bnd2++;
67 }
68 }
69 bnd1++;
70 }
71
72 for (mci = marker.begin(); mci != marker.end(); mci++) {
73 charges[*mci] = (formal / marker.size());
74 }
75 }
76 }
77 /*
78 for (aix = 0; aix < natms; aix++) {
79 std::cout << "In splitter: " << " charges:" << charges[aix] << "\n";
80 }*/
81 }
82 } // end of namespace Gasteiger
83
84 namespace RDKit {
computeGasteigerCharges(const ROMol * mol,int nIter,bool throwOnParamFailure)85 void computeGasteigerCharges(const ROMol *mol, int nIter,
86 bool throwOnParamFailure) {
87 PRECONDITION(mol, "bad molecule");
88 computeGasteigerCharges(*mol, nIter, throwOnParamFailure);
89 }
computeGasteigerCharges(const ROMol & mol,int nIter,bool throwOnParamFailure)90 void computeGasteigerCharges(const ROMol &mol, int nIter,
91 bool throwOnParamFailure) {
92 std::vector<double> chgs(mol.getNumAtoms());
93 computeGasteigerCharges(mol, chgs, nIter, throwOnParamFailure);
94 }
95
96 /*! \brief compute the Gasteiger partial charges and return a new molecule with
97 *the charges set
98 *
99 * Ref : J.Gasteiger, M. Marsili, "Iterative Equalization of Orbital
100 * Electronegativity
101 * A Rapid Access to Atomic Charges", Tetrahedron Vol 36 p3219 1980
102 */
computeGasteigerCharges(const ROMol & mol,std::vector<double> & charges,int nIter,bool throwOnParamFailure)103 void computeGasteigerCharges(const ROMol &mol, std::vector<double> &charges,
104 int nIter, bool throwOnParamFailure) {
105 PRECONDITION(charges.size() >= mol.getNumAtoms(), "bad array size");
106
107 PeriodicTable *table = PeriodicTable::getTable();
108 const GasteigerParams *params = GasteigerParams::getParams();
109
110 double damp = DAMP;
111 int natms = mol.getNumAtoms();
112 // space for parameters for each atom in the molecule
113 std::vector<DOUBLE_VECT> atmPs;
114 atmPs.reserve(natms);
115
116 std::fill(charges.begin(), charges.end(), 0.0);
117
118 DOUBLE_VECT
119 hChrg; // total charge on the implicit hydrogen on each heavy atom
120 hChrg.resize(natms, 0.0);
121
122 DOUBLE_VECT ionX;
123 ionX.resize(natms, 0.0);
124
125 DOUBLE_VECT energ;
126 energ.resize(natms, 0.0);
127
128 ROMol::ADJ_ITER nbrIdx, endIdx;
129
130 // deal with the conjugated system - distribute the formal charges on atoms of
131 // same type in each
132 // conjugated system
133 Gasteiger::splitChargeConjugated(mol, charges);
134
135 // now read in the parameters
136 ROMol::ConstAtomIterator ai;
137
138 for (ai = mol.beginAtoms(); ai != mol.endAtoms(); ai++) {
139 std::string elem = table->getElementSymbol((*ai)->getAtomicNum());
140 std::string mode;
141
142 switch ((*ai)->getHybridization()) {
143 case Atom::SP3:
144 mode = "sp3";
145 break;
146 case Atom::SP2:
147 mode = "sp2";
148 break;
149 case Atom::SP:
150 mode = "sp";
151 break;
152 default:
153 if ((*ai)->getAtomicNum() == 1) {
154 // if it is hydrogen
155 mode = "*";
156 } else if ((*ai)->getAtomicNum() == 16) {
157 // we have a sulfur atom with no hybridization information
158 // check how many oxygens we have on the sulfur
159 boost::tie(nbrIdx, endIdx) = mol.getAtomNeighbors(*ai);
160 int no = 0;
161 while (nbrIdx != endIdx) {
162 if (mol.getAtomWithIdx(*nbrIdx)->getAtomicNum() == 8) {
163 no++;
164 }
165 nbrIdx++;
166 }
167 if (no == 2) {
168 mode = "so2";
169 } else if (no == 1) {
170 mode = "so";
171 } else {
172 // some other sulfur state. Default to sp3
173 mode = "sp3";
174 }
175 }
176 }
177
178 // if we get a unknown mode or element type the
179 // following will will throw an exception
180 atmPs.push_back(params->getParams(elem, mode, throwOnParamFailure));
181
182 // set ionX parameters
183 // if Hydrogen treat differently
184 int idx = (*ai)->getIdx();
185 if ((*ai)->getAtomicNum() == 1) {
186 ionX[idx] = IONXH;
187 } else {
188 ionX[idx] = atmPs[idx][0] + atmPs[idx][1] + atmPs[idx][2];
189 }
190 }
191
192 // do the iteration here
193 int itx, aix, sgn, niHs;
194 double enr, dq, dx, qHs, dqH;
195 // parameters for hydrogen atoms (for case where the hydrogen are not in the
196 // graph (implicit hydrogens)
197 DOUBLE_VECT hParams;
198 hParams = params->getParams("H", "*", throwOnParamFailure);
199
200 /*
201 int itmp;
202 for (itmp = 0; itmp < 5; itmp++) {
203 std::cout << " aq:" << charges[itmp] << "\n";
204 }*/
205
206 for (itx = 0; itx < nIter; itx++) {
207 for (aix = 0; aix < natms; aix++) {
208 // enr = p0 + charge*(p1 + p2*charge)
209 enr = atmPs[aix][0] +
210 charges[aix] * (atmPs[aix][1] + atmPs[aix][2] * charges[aix]);
211 energ[aix] = enr;
212 }
213
214 for (aix = 0; aix < natms; aix++) {
215 dq = 0.0;
216 boost::tie(nbrIdx, endIdx) =
217 mol.getAtomNeighbors(mol.getAtomWithIdx(aix));
218 while (nbrIdx != endIdx) {
219 dx = energ[*nbrIdx] - energ[aix];
220 if (dx < 0.0) {
221 sgn = 0;
222 } else {
223 sgn = 1;
224 }
225 dq += dx / ((sgn * (ionX[aix] - ionX[*nbrIdx])) + ionX[*nbrIdx]);
226 nbrIdx++;
227 }
228 // now loop over the implicit hydrogens and get their contributions
229 // since hydrogens don't connect to anything else, update their charges at
230 // the same time
231 niHs = mol.getAtomWithIdx(aix)->getTotalNumHs();
232 if (niHs > 0) {
233 qHs = hChrg[aix] / niHs;
234 enr = hParams[0] + qHs * (hParams[1] + hParams[2] * qHs);
235 dx = enr - energ[aix];
236 if (dx < 0.0) {
237 sgn = 0;
238 } else {
239 sgn = 1;
240 }
241
242 dqH = dx / ((sgn * (ionX[aix] - IONXH)) + IONXH);
243
244 dq += (niHs * dqH);
245
246 // adjust the charges on the hydrogens simultaneously (possible because
247 // each of the
248 // hydrogens have no other neighbors)
249 hChrg[aix] -= (niHs * dqH * damp);
250 }
251 charges[aix] += (damp * dq);
252 }
253
254 damp *= DAMP_SCALE;
255 }
256
257 for (aix = 0; aix < natms; aix++) {
258 mol.getAtomWithIdx(aix)->setProp(common_properties::_GasteigerCharge,
259 charges[aix], true);
260 // set the implicit hydrogen charges
261 mol.getAtomWithIdx(aix)->setProp(common_properties::_GasteigerHCharge,
262 hChrg[aix], true);
263 }
264 }
265 }
266