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