src/openbabel/molchrg.cpp

/**********************************************************************
molchrg.cpp - Assign Gasteiger partial charges.
 
Copyright (C) 1998-2001 by OpenEye Scientific Software, Inc.
Some portions Copyright (C) 2001-2005 by Geoffrey R. Hutchison
 
This file is part of the Open Babel project.
For more information, see <http://openbabel.sourceforge.net/>
 
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation version 2 of the License.
 
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
***********************************************************************/

#include "mol.hpp"
#include "molchrg.hpp"

using namespace std;
namespace OpenBabel
{

/*! \class OBGastChrg
    \brief Assigns Gasteiger partial charges

The OBGastChrg class is responsible for the assignment of partial
charges to a molecule according to the Gasteiger charge model (sigma). When
the partial charge of an atom is requested and partial charges do not
yet exist for the molecule the OBGastChrg class will automatically be
called to assign charges for the molecule. If charges have been read
in for a molecule the following code shows how to force the
recalculation of partial charges:
\code
OBMol mol;

mol.UnsetPartialChargesPerceived();
FOR_ATOMS_IN_MOL(atom, mol)
{
 cout << "atom number = " << atom->GetIdx();
 cout << " charge = " << atom->GetPartialCharge() << endl;
}
\endcode
Formal charges are used as seed values of the initial charge of atoms,
and the partial charge is propagated to neighboring atoms. For
example, quaternary amines would have a +1 charge, and the effect of
the positive charge would be felt by neighboring atoms according to
the Gasteiger model (sigma).

For more information, see:
J. Gasteiger & M. Marsili, "A New Model for Calculating Atomic Charges in Molecules" Tetrahedron Lett., (1978) 3181-3184.
*/

bool OBGastChrg::AssignPartialCharges(OBMol &mol)
{
    //InitialPartialCharges(mol);
    obErrorLog.ThrowError(__FUNCTION__,
                          "Ran OpenBabel::AssignPartialCharges", obAuditMsg);

    OBAtom *atom;
    vector<OBNodeBase*>::iterator i;

    GSVResize(mol.NumAtoms()+1);
    double a,b,c;
    for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i))
    {
        if (!GasteigerSigmaChi(atom,a,b,c))
            return(false);
        _gsv[atom->GetIdx()]->SetValues(a,b,c,atom->GetPartialCharge());
    }

    double alpha,charge,denom;
    unsigned j;
    int iter;
    OBBond *bond;
    OBAtom *src,*dst;
    vector<OBEdgeBase*>::iterator k;
    alpha = 1.0;
    for(iter = 0;iter < OB_GASTEIGER_ITERS;iter++)
    {
        alpha *= OB_GASTEIGER_DAMP;

        for( j=1;j < _gsv.size();j++)
        {
            charge = _gsv[j]->q;
            _gsv[j]->chi = (_gsv[j]->c*charge+_gsv[j]->b)*charge+_gsv[j]->a;
        }

        for (bond = mol.BeginBond(k);bond;bond = mol.NextBond(k))
        {
            src = bond->GetBeginAtom();
            dst = bond->GetEndAtom();

            if (_gsv[src->GetIdx()]->chi >= _gsv[dst->GetIdx()]->chi)
            {
                if (dst->IsHydrogen())
                    denom = double(OB_GASTEIGER_DENOM);
                else
                    denom = _gsv[dst->GetIdx()]->denom;
            }
            else
            {
                if (src->IsHydrogen())
                    denom = double(OB_GASTEIGER_DENOM);
                else
                    denom = _gsv[src->GetIdx()]->denom;
            }

            charge = (_gsv[src->GetIdx()]->chi - _gsv[dst->GetIdx()]->chi)/denom;
            _gsv[src->GetIdx()]->q -= alpha*charge;
            _gsv[dst->GetIdx()]->q += alpha*charge;
        }
    }

    for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i))
        atom->SetPartialCharge(_gsv[atom->GetIdx()]->q);

    return(true);
}

void OBGastChrg::InitialPartialCharges(OBMol &mol)
{
    OBAtom *atom;
    vector<OBNodeBase*>::iterator i;

    for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i))
    {
        if (atom->IsCarboxylOxygen())
            atom->SetPartialCharge(-0.500);
        else if (atom->IsPhosphateOxygen() &&
                 atom->GetHvyValence() == 1)
            atom->SetPartialCharge(-0.666);
        else if (atom->IsSulfateOxygen())
            atom->SetPartialCharge(-0.500);
        else
            atom->SetPartialCharge((double)atom->GetFormalCharge());
    }
}

bool OBGastChrg::GasteigerSigmaChi(OBAtom *atom,double &a,double &b,double &c )
{
    int count;
    double val[3] = {0.0,0.0,0.0};

    switch(atom->GetAtomicNum())
    {
    case 1: //H
        val[0] = 0.37;
        val[1] = 7.17;
        val[2] = 12.85;
        break;
    case 6: //C
        if (atom->GetHyb() == 3)
        {
            val[0] = 0.68;
            val[1] = 7.98;
            val[2] = 19.04;
        }
        if (atom->GetHyb() == 2)
        {
            val[0] = 0.98;
            val[1] = 8.79;
            val[2] = 19.62;
        }
        if (atom->GetHyb() == 1)
        {
            val[0] = 1.67;
            val[1] = 10.39;
            val[2] = 20.57;
        }
        break;
    case 7: //N
        if (atom->GetHyb() == 3)
        {
            if (atom->GetValence() == 4 || atom->GetFormalCharge())
            {
                val[0] = 0.0;
                val[1] = 0.0;
                val[2] = 23.72;
            }
            else
            {
                val[0] = 2.08;
                val[1] = 11.54;
                val[2] = 23.72;
            }
        }

        if (atom->GetHyb() == 2)
        {
            if (EQ(atom->GetType(),"Npl") || EQ(atom->GetType(),"Nam"))
            {
                val[0] = 2.46;
                val[1] = 12.32;
                val[2] = 24.86;
            }
            else
            {
                val[0] = 2.57;
                val[1] = 12.87;
                val[2] = 24.87;
            }
        }

        if (atom->GetHyb() == 1)
        {
            val[0] = 3.71;
            val[1] = 15.68;
            val[2] = 27.11;
        }
        break;
    case 8: //O
        if (atom->GetHyb() == 3)
        {
            val[0] = 2.65;
            val[1] = 14.18;
            val[2] = 28.49;
        }
        if (atom->GetHyb() == 2)
        {
            val[0] = 3.75;
            val[1] = 17.07;
            val[2] = 31.33;
        }
        break;
    case 9: //F
        val[0] = 3.12;
        val[1] = 14.66;
        val[2] = 30.82;
        break;
    case 15: //P
        val[0] = 1.62;
        val[1] = 8.90;
        val[2] = 18.10;
        break;
    case 16: //S
        count = atom->CountFreeOxygens();
        if (count == 0 || count == 1)
        {
            val[0] = 2.39;
            val[1] = 10.14;
            val[2] = 20.65;
        }
        if (count > 1)
        {
            val[0] = 2.39;
            val[1] = 12.00;
            val[2] = 24.00;
        }
        /*S2? if (count == 0) {val[0] = 2.72;val[1] = 10.88;val[2] = 21.69;}*/
        break;
    case 17: //Cl
        val[0] = 2.66;
        val[1] = 11.00;
        val[2] = 22.04;
        break;
    case 35: //Br
        val[0] = 2.77;
        val[1] = 10.08;
        val[2] = 19.71;
        break;
    case 53: //I
        val[0] = 2.90;
        val[1] = 9.90;
        val[2] = 18.82;
        break;
    case 13: //Al
        val[0] = 1.06;
        val[1] = 5.47;
        val[2] = 11.65;
        break;
    }

    if (!IsNearZero(val[2]))
    {
        a = val[1];
        b = (val[2]-val[0])/2;
        c = (val[2]+val[0])/2 - val[1];
    }
    else
        return(false);

    return(true);
}

void OBGastChrg::GSVResize(int size)
{
    vector <GasteigerState*>::iterator i;
    for (i = _gsv.begin();i != _gsv.end();i++)
        delete *i;
    _gsv.clear();

    for (int j = 0;j < size;j++)
        _gsv.push_back(new GasteigerState);
}

OBGastChrg::~OBGastChrg()
{
    vector <GasteigerState*>::iterator i;
    for (i = _gsv.begin();i != _gsv.end();i++)
        delete *i;
}

GasteigerState::GasteigerState()
{
    a = 0.0;
    b = 0.0;
    c = 0.0;
    denom = 0.0;
    chi = 0.0;
    q = 0.0;
}

} // end namespace OpenBabel

//! \file molchrg.cpp
//! \brief Assign Gasteiger partial charges.
Revision:	741
Committed:	Wed Nov 16 19:42:11 2005 UTC (19 years, 5 months ago) by tim
File size:	8405 byte(s)
Log Message:	adding openbabel
#	User	Rev	Content
1	tim	741	/**********************************************************************
2			molchrg.cpp - Assign Gasteiger partial charges.
3
4			Copyright (C) 1998-2001 by OpenEye Scientific Software, Inc.
5			Some portions Copyright (C) 2001-2005 by Geoffrey R. Hutchison
6
7			This file is part of the Open Babel project.
8			For more information, see <http://openbabel.sourceforge.net/>
9
10			This program is free software; you can redistribute it and/or modify
11			it under the terms of the GNU General Public License as published by
12			the Free Software Foundation version 2 of the License.
13
14			This program is distributed in the hope that it will be useful,
15			but WITHOUT ANY WARRANTY; without even the implied warranty of
16			MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17			GNU General Public License for more details.
18			***********************************************************************/
19
20			#include "mol.hpp"
21			#include "molchrg.hpp"
22
23			using namespace std;
24			namespace OpenBabel
25			{
26
27			/*! \class OBGastChrg
28			\brief Assigns Gasteiger partial charges
29
30			The OBGastChrg class is responsible for the assignment of partial
31			charges to a molecule according to the Gasteiger charge model (sigma). When
32			the partial charge of an atom is requested and partial charges do not
33			yet exist for the molecule the OBGastChrg class will automatically be
34			called to assign charges for the molecule. If charges have been read
35			in for a molecule the following code shows how to force the
36			recalculation of partial charges:
37			\code
38			OBMol mol;
39
40			mol.UnsetPartialChargesPerceived();
41			FOR_ATOMS_IN_MOL(atom, mol)
42			{
43			cout << "atom number = " << atom->GetIdx();
44			cout << " charge = " << atom->GetPartialCharge() << endl;
45			}
46			\endcode
47			Formal charges are used as seed values of the initial charge of atoms,
48			and the partial charge is propagated to neighboring atoms. For
49			example, quaternary amines would have a +1 charge, and the effect of
50			the positive charge would be felt by neighboring atoms according to
51			the Gasteiger model (sigma).
52
53			For more information, see:
54			J. Gasteiger & M. Marsili, "A New Model for Calculating Atomic Charges in Molecules" Tetrahedron Lett., (1978) 3181-3184.
55			*/
56
57			bool OBGastChrg::AssignPartialCharges(OBMol &mol)
58			{
59			//InitialPartialCharges(mol);
60			obErrorLog.ThrowError(__FUNCTION__,
61			"Ran OpenBabel::AssignPartialCharges", obAuditMsg);
62
63			OBAtom *atom;
64			vector<OBNodeBase*>::iterator i;
65
66			GSVResize(mol.NumAtoms()+1);
67			double a,b,c;
68			for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i))
69			{
70			if (!GasteigerSigmaChi(atom,a,b,c))
71			return(false);
72			_gsv[atom->GetIdx()]->SetValues(a,b,c,atom->GetPartialCharge());
73			}
74
75			double alpha,charge,denom;
76			unsigned j;
77			int iter;
78			OBBond *bond;
79			OBAtom src,dst;
80			vector<OBEdgeBase*>::iterator k;
81			alpha = 1.0;
82			for(iter = 0;iter < OB_GASTEIGER_ITERS;iter++)
83			{
84			alpha *= OB_GASTEIGER_DAMP;
85
86			for( j=1;j < _gsv.size();j++)
87			{
88			charge = _gsv[j]->q;
89			_gsv[j]->chi = (_gsv[j]->ccharge+_gsv[j]->b)charge+_gsv[j]->a;
90			}
91
92			for (bond = mol.BeginBond(k);bond;bond = mol.NextBond(k))
93			{
94			src = bond->GetBeginAtom();
95			dst = bond->GetEndAtom();
96
97			if (_gsv[src->GetIdx()]->chi >= _gsv[dst->GetIdx()]->chi)
98			{
99			if (dst->IsHydrogen())
100			denom = double(OB_GASTEIGER_DENOM);
101			else
102			denom = _gsv[dst->GetIdx()]->denom;
103			}
104			else
105			{
106			if (src->IsHydrogen())
107			denom = double(OB_GASTEIGER_DENOM);
108			else
109			denom = _gsv[src->GetIdx()]->denom;
110			}
111
112			charge = (_gsv[src->GetIdx()]->chi - _gsv[dst->GetIdx()]->chi)/denom;
113			_gsv[src->GetIdx()]->q -= alpha*charge;
114			_gsv[dst->GetIdx()]->q += alpha*charge;
115			}
116			}
117
118			for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i))
119			atom->SetPartialCharge(_gsv[atom->GetIdx()]->q);
120
121			return(true);
122			}
123
124			void OBGastChrg::InitialPartialCharges(OBMol &mol)
125			{
126			OBAtom *atom;
127			vector<OBNodeBase*>::iterator i;
128
129			for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i))
130			{
131			if (atom->IsCarboxylOxygen())
132			atom->SetPartialCharge(-0.500);
133			else if (atom->IsPhosphateOxygen() &&
134			atom->GetHvyValence() == 1)
135			atom->SetPartialCharge(-0.666);
136			else if (atom->IsSulfateOxygen())
137			atom->SetPartialCharge(-0.500);
138			else
139			atom->SetPartialCharge((double)atom->GetFormalCharge());
140			}
141			}
142
143			bool OBGastChrg::GasteigerSigmaChi(OBAtom *atom,double &a,double &b,double &c )
144			{
145			int count;
146			double val[3] = {0.0,0.0,0.0};
147
148			switch(atom->GetAtomicNum())
149			{
150			case 1: //H
151			val[0] = 0.37;
152			val[1] = 7.17;
153			val[2] = 12.85;
154			break;
155			case 6: //C
156			if (atom->GetHyb() == 3)
157			{
158			val[0] = 0.68;
159			val[1] = 7.98;
160			val[2] = 19.04;
161			}
162			if (atom->GetHyb() == 2)
163			{
164			val[0] = 0.98;
165			val[1] = 8.79;
166			val[2] = 19.62;
167			}
168			if (atom->GetHyb() == 1)
169			{
170			val[0] = 1.67;
171			val[1] = 10.39;
172			val[2] = 20.57;
173			}
174			break;
175			case 7: //N
176			if (atom->GetHyb() == 3)
177			{
178			if (atom->GetValence() == 4 \|\| atom->GetFormalCharge())
179			{
180			val[0] = 0.0;
181			val[1] = 0.0;
182			val[2] = 23.72;
183			}
184			else
185			{
186			val[0] = 2.08;
187			val[1] = 11.54;
188			val[2] = 23.72;
189			}
190			}
191
192			if (atom->GetHyb() == 2)
193			{
194			if (EQ(atom->GetType(),"Npl") \|\| EQ(atom->GetType(),"Nam"))
195			{
196			val[0] = 2.46;
197			val[1] = 12.32;
198			val[2] = 24.86;
199			}
200			else
201			{
202			val[0] = 2.57;
203			val[1] = 12.87;
204			val[2] = 24.87;
205			}
206			}
207
208			if (atom->GetHyb() == 1)
209			{
210			val[0] = 3.71;
211			val[1] = 15.68;
212			val[2] = 27.11;
213			}
214			break;
215			case 8: //O
216			if (atom->GetHyb() == 3)
217			{
218			val[0] = 2.65;
219			val[1] = 14.18;
220			val[2] = 28.49;
221			}
222			if (atom->GetHyb() == 2)
223			{
224			val[0] = 3.75;
225			val[1] = 17.07;
226			val[2] = 31.33;
227			}
228			break;
229			case 9: //F
230			val[0] = 3.12;
231			val[1] = 14.66;
232			val[2] = 30.82;
233			break;
234			case 15: //P
235			val[0] = 1.62;
236			val[1] = 8.90;
237			val[2] = 18.10;
238			break;
239			case 16: //S
240			count = atom->CountFreeOxygens();
241			if (count == 0 \|\| count == 1)
242			{
243			val[0] = 2.39;
244			val[1] = 10.14;
245			val[2] = 20.65;
246			}
247			if (count > 1)
248			{
249			val[0] = 2.39;
250			val[1] = 12.00;
251			val[2] = 24.00;
252			}
253			/S2? if (count == 0) {val[0] = 2.72;val[1] = 10.88;val[2] = 21.69;}/
254			break;
255			case 17: //Cl
256			val[0] = 2.66;
257			val[1] = 11.00;
258			val[2] = 22.04;
259			break;
260			case 35: //Br
261			val[0] = 2.77;
262			val[1] = 10.08;
263			val[2] = 19.71;
264			break;
265			case 53: //I
266			val[0] = 2.90;
267			val[1] = 9.90;
268			val[2] = 18.82;
269			break;
270			case 13: //Al
271			val[0] = 1.06;
272			val[1] = 5.47;
273			val[2] = 11.65;
274			break;
275			}
276
277			if (!IsNearZero(val[2]))
278			{
279			a = val[1];
280			b = (val[2]-val[0])/2;
281			c = (val[2]+val[0])/2 - val[1];
282			}
283			else
284			return(false);
285
286			return(true);
287			}
288
289			void OBGastChrg::GSVResize(int size)
290			{
291			vector <GasteigerState*>::iterator i;
292			for (i = _gsv.begin();i != _gsv.end();i++)
293			delete *i;
294			_gsv.clear();
295
296			for (int j = 0;j < size;j++)
297			_gsv.push_back(new GasteigerState);
298			}
299
300			OBGastChrg::~OBGastChrg()
301			{
302			vector <GasteigerState*>::iterator i;
303			for (i = _gsv.begin();i != _gsv.end();i++)
304			delete *i;
305			}
306
307			GasteigerState::GasteigerState()
308			{
309			a = 0.0;
310			b = 0.0;
311			c = 0.0;
312			denom = 0.0;
313			chi = 0.0;
314			q = 0.0;
315			}
316
317			} // end namespace OpenBabel
318
319			//! \file molchrg.cpp
320			//! \brief Assign Gasteiger partial charges.