src/openbabel/phmodel.cpp

/**********************************************************************
phmodel.cpp - Read pH rules and assign 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 "phmodel.hpp"
#include "phmodeldata.hpp"

#ifdef WIN32
#pragma warning (disable : 4786)
#endif

using namespace std;

namespace OpenBabel
{

  // Global OBPhModel for assigning formal charges and hydrogen addition rules
OBPhModel phmodel;
extern OBAtomTyper atomtyper;

OBPhModel::OBPhModel()
{
    _init = false;
    STR_DEFINE(_dir, FRC_PATH);
    _envvar = "FORCE_PARAM_PATH";
    _filename = "phmodel.txt";
    _subdir = "data";
    _dataptr = PhModelData;
}

OBPhModel::~OBPhModel()
{
    vector<OBChemTsfm*>::iterator k;
    for (k = _vtsfm.begin();k != _vtsfm.end();k++)
        delete *k;

    vector<pair<OBSmartsPattern*,vector<double> > >::iterator m;
    for (m = _vschrg.begin();m != _vschrg.end();m++)
        delete m->first;
}

void OBPhModel::ParseLine(const char *buffer)
{
    vector<string> vs;
    OBSmartsPattern *sp;

    if (buffer[0] == '#')
        return;

    if (EQn(buffer,"TRANSFORM",7))
    {
        tokenize(vs,buffer);
        if (vs.empty() || vs.size() < 4)
          {
            obErrorLog.ThrowError(__FUNCTION__, " Could not parse line in phmodel table from phmodel.txt", obInfo);
            return;
          }

        OBChemTsfm *tsfm = new OBChemTsfm;
        if (!tsfm->Init(vs[1],vs[3]))
        {
            delete tsfm;
            tsfm = NULL;
            obErrorLog.ThrowError(__FUNCTION__, " Could not parse line in phmodel table from phmodel.txt", obInfo);
            return;
        }

        _vtsfm.push_back(tsfm);
    }
    else if (EQn(buffer,"SEEDCHARGE",10))
    {
        tokenize(vs,buffer);
        if (vs.empty() || vs.size() < 2)
          {
            obErrorLog.ThrowError(__FUNCTION__, " Could not parse line in phmodel table from phmodel.txt", obInfo);
            return;
          }

        sp = new OBSmartsPattern;
        if (!sp->Init(vs[1]) || (vs.size()-2) != sp->NumAtoms())
        {
            delete sp;
            sp = NULL;
            obErrorLog.ThrowError(__FUNCTION__, " Could not parse line in phmodel table from phmodel.txt", obInfo);
            return;
        }

        vector<double> vf;
        vector<string>::iterator i;
        for (i = vs.begin()+2;i != vs.end();i++)
            vf.push_back(atof((char*)i->c_str()));

        _vschrg.push_back(pair<OBSmartsPattern*,vector<double> > (sp,vf));
    }
}

void OBPhModel::AssignSeedPartialCharge(OBMol &mol)
{
    if (!_init)
        Init();

    mol.SetPartialChargesPerceived();
    if (!mol.AutomaticPartialCharge())
        return;

    vector<pair<OBSmartsPattern*,vector<double> > >::iterator i;
    for (i = _vschrg.begin();i != _vschrg.end();i++)
        if (i->first->Match(mol))
        {
            _mlist = i->first->GetUMapList();
            unsigned int k;
            vector<vector<int> >::iterator j;

            for (j = _mlist.begin();j != _mlist.end();j++)
                for (k = 0;k < j->size();k++)
                    mol.GetAtom((*j)[k])->SetPartialCharge(i->second[k]);
        }
}

void OBPhModel::CorrectForPH(OBMol &mol)
{
    if (!_init)
        Init();
    if (mol.IsCorrectedForPH())
        return;
    if (!mol.AutomaticFormalCharge())
        return;

    mol.SetCorrectedForPH();

    obErrorLog.ThrowError(__FUNCTION__,
                          "Ran OpenBabel::CorrectForPH", obAuditMsg);

    OBAtom *atom;
    vector<OBNodeBase*>::iterator j;
    for (atom = mol.BeginAtom(j);atom;atom = mol.NextAtom(j))
        atom->SetFormalCharge(0);

    vector<OBChemTsfm*>::iterator i;
    for (i = _vtsfm.begin();i != _vtsfm.end();i++)
        (*i)->Apply(mol);

    atomtyper.CorrectAromaticNitrogens(mol);
}

// Portions of this documentation adapted from the JOELib docs, written by
// Joerg Wegner
/** \class OBChemTsfm
    \brief SMARTS based structural modification (chemical transformation)
 
    Transformation of chemical structures can be used for pH value correction
    (i.e. via OBPhModel and OBMol::CorrectForPH()). The OBChemTsfm class
    defines SMARTS based TRANSFORM patterns to delete atoms, change atom types,
    atom formal charges, and bond types.

    For storing and converting chemical reaction files, use the OBReaction class.
 **/
bool OBChemTsfm::Init(string &bgn,string &end)
{
    if (!_bgn.Init(bgn))
        return(false);
    if (!end.empty())
        if (!_end.Init(end))
            return(false);

    //find atoms to be deleted
    unsigned int i,j;
    int vb;
    bool found;
    for (i = 0;i < _bgn.NumAtoms();i++)
        if ((vb = _bgn.GetVectorBinding(i)))
        {
            found = false;
            for (j = 0;j < _end.NumAtoms();j++)
                if (vb == _end.GetVectorBinding(j))
                {
                    found = true;
                    break;
                }

            if (!found)
              _vadel.push_back(i);
        }

    //find elements to be changed
    int ele;
    for (i = 0;i < _bgn.NumAtoms();i++)
        if ((vb = _bgn.GetVectorBinding(i)) != 0)
        {
            ele = _bgn.GetAtomicNum(i);
            for (j = 0;j < _end.NumAtoms();j++)
                if (vb == _end.GetVectorBinding(j))
                    if (ele != _end.GetAtomicNum(j))
                    {
                        _vele.push_back(pair<int,int> (i,_end.GetAtomicNum(j)));
                        break;
                    }
        }

    //find charges to modify
    int chrg;
    for (i = 0;i < _bgn.NumAtoms();i++)
        if ((vb = _bgn.GetVectorBinding(i)))
        {
            chrg = _bgn.GetCharge(i);
            for (j = 0;j < _end.NumAtoms();j++)
                if (vb == _end.GetVectorBinding(j))
                    if (chrg != _end.GetCharge(j))
                        _vchrg.push_back(pair<int,int> (i,_end.GetCharge(j)));
        }

    //find bonds to be modified
    //find bonds to be modified
    int bsrc,bdst,bord,bvb1,bvb2;
    int esrc,edst,eord,evb1,evb2;

    for (i = 0;i < _bgn.NumBonds();i++)
    {
        _bgn.GetBond(bsrc,bdst,bord,i);
        bvb1 = _bgn.GetVectorBinding(bsrc);
        bvb2 = _bgn.GetVectorBinding(bdst);
        if (!bvb1 || !bvb2)
            continue;

        for (j = 0;j < _end.NumBonds();j++)
        {
            _end.GetBond(esrc,edst,eord,j);
            evb1 = _end.GetVectorBinding(esrc);
            evb2 = _end.GetVectorBinding(edst);
            if ((bvb1 == evb1 && bvb2 == evb2) || (bvb1 == evb2 && bvb2 == evb1))
            {
                if (bord == eord)
                    break; //nothing to modify if bond orders identical
                _vbond.push_back(pair<pair<int,int>,int> (pair<int,int> (bsrc,bdst),eord));
                break;
            }
        }
    }

    //make sure there is some kind of transform to do here
    if (_vadel.empty() && _vchrg.empty() && _vbond.empty())
        return(false);

    return(true);
}

bool OBChemTsfm::Apply(OBMol &mol)
{
    if (!_bgn.Match(mol))
        return(false);

    vector<vector<int> > mlist = _bgn.GetUMapList();

    obErrorLog.ThrowError(__FUNCTION__,
                          "Ran OpenBabel::OBChemTransform", obAuditMsg);

    if (!_vchrg.empty()) //modify charges
    {
        vector<vector<int> >::iterator i;
        vector<pair<int,int> >::iterator j;

        for (i = mlist.begin();i != mlist.end();i++)
            for (j = _vchrg.begin();j != _vchrg.end();j++)
                if (j->first < (signed)i->size()) //goof proofing
                    mol.GetAtom((*i)[j->first])->SetFormalCharge(j->second);

        mol.UnsetImplicitValencePerceived();
    }

    if (!_vbond.empty()) //modify bond orders
    {
        OBBond *bond;
        vector<vector<int> >::iterator i;
        vector<pair<pair<int,int>,int> >::iterator j;
        for (i = mlist.begin();i != mlist.end();i++)
            for (j = _vbond.begin();j != _vbond.end();j++)
            {
                bond = mol.GetBond((*i)[j->first.first],(*i)[j->first.second]);
                if (!bond)
                {
                  obErrorLog.ThrowError(__FUNCTION__, "unable to find bond", obDebug);
                    continue;
                }

                bond->SetBO(j->second);
            }
    }

    if (!_vadel.empty() || !_vele.empty()) //delete atoms and change elements
    {
        vector<int>::iterator j;
        vector<vector<int> >::iterator i;

        if (!_vele.empty())
        {
            vector<pair<int,int> >::iterator k;
            for (i = mlist.begin();i != mlist.end();i++)
                for (k = _vele.begin();k != _vele.end();k++)
                    mol.GetAtom((*i)[k->first])->SetAtomicNum(k->second);
        }

        //make sure same atom isn't deleted twice
        vector<bool> vda;
        vector<OBNodeBase*> vdel;
        vda.resize(mol.NumAtoms()+1,false);
        for (i = mlist.begin();i != mlist.end();i++)
            for (j = _vadel.begin();j != _vadel.end();j++)
                if (!vda[(*i)[*j]])
                {
                    vda[(*i)[*j]] = true;
                    vdel.push_back(mol.GetAtom((*i)[*j]));
                }

        vector<OBNodeBase*>::iterator k;
        for (k = vdel.begin();k != vdel.end();k++)
            mol.DeleteAtom((OBAtom*)*k);
    }

    return(true);
}

} //namespace OpenBabel

//! \file phmodel.cpp
//! \brief Read pH rules and assign charges.
Revision:	751
Committed:	Thu Nov 17 20:38:45 2005 UTC (19 years, 11 months ago) by gezelter
File size:	10054 byte(s)
Log Message:	Unifying config.h stuff and making sure the OpenBabel codes can find our default (and environment variable) Force Field directories.
#	User	Rev	Content
1	tim	741	/**********************************************************************
2			phmodel.cpp - Read pH rules and assign 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 "phmodel.hpp"
22			#include "phmodeldata.hpp"
23
24			#ifdef WIN32
25			#pragma warning (disable : 4786)
26			#endif
27
28			using namespace std;
29
30			namespace OpenBabel
31			{
32
33			// Global OBPhModel for assigning formal charges and hydrogen addition rules
34			OBPhModel phmodel;
35			extern OBAtomTyper atomtyper;
36
37			OBPhModel::OBPhModel()
38			{
39			_init = false;
40	gezelter	751	STR_DEFINE(_dir, FRC_PATH);
41			_envvar = "FORCE_PARAM_PATH";
42	tim	741	_filename = "phmodel.txt";
43			_subdir = "data";
44			_dataptr = PhModelData;
45			}
46
47			OBPhModel::~OBPhModel()
48			{
49			vector<OBChemTsfm*>::iterator k;
50			for (k = _vtsfm.begin();k != _vtsfm.end();k++)
51			delete *k;
52
53			vector<pair<OBSmartsPattern*,vector<double> > >::iterator m;
54			for (m = _vschrg.begin();m != _vschrg.end();m++)
55			delete m->first;
56			}
57
58			void OBPhModel::ParseLine(const char *buffer)
59			{
60			vector<string> vs;
61			OBSmartsPattern *sp;
62
63			if (buffer[0] == '#')
64			return;
65
66			if (EQn(buffer,"TRANSFORM",7))
67			{
68			tokenize(vs,buffer);
69			if (vs.empty() \|\| vs.size() < 4)
70			{
71			obErrorLog.ThrowError(__FUNCTION__, " Could not parse line in phmodel table from phmodel.txt", obInfo);
72			return;
73			}
74
75			OBChemTsfm *tsfm = new OBChemTsfm;
76			if (!tsfm->Init(vs[1],vs[3]))
77			{
78			delete tsfm;
79			tsfm = NULL;
80			obErrorLog.ThrowError(__FUNCTION__, " Could not parse line in phmodel table from phmodel.txt", obInfo);
81			return;
82			}
83
84			_vtsfm.push_back(tsfm);
85			}
86			else if (EQn(buffer,"SEEDCHARGE",10))
87			{
88			tokenize(vs,buffer);
89			if (vs.empty() \|\| vs.size() < 2)
90			{
91			obErrorLog.ThrowError(__FUNCTION__, " Could not parse line in phmodel table from phmodel.txt", obInfo);
92			return;
93			}
94
95			sp = new OBSmartsPattern;
96			if (!sp->Init(vs[1]) \|\| (vs.size()-2) != sp->NumAtoms())
97			{
98			delete sp;
99			sp = NULL;
100			obErrorLog.ThrowError(__FUNCTION__, " Could not parse line in phmodel table from phmodel.txt", obInfo);
101			return;
102			}
103
104			vector<double> vf;
105			vector<string>::iterator i;
106			for (i = vs.begin()+2;i != vs.end();i++)
107			vf.push_back(atof((char*)i->c_str()));
108
109			_vschrg.push_back(pair<OBSmartsPattern*,vector<double> > (sp,vf));
110			}
111			}
112
113			void OBPhModel::AssignSeedPartialCharge(OBMol &mol)
114			{
115			if (!_init)
116			Init();
117
118			mol.SetPartialChargesPerceived();
119			if (!mol.AutomaticPartialCharge())
120			return;
121
122			vector<pair<OBSmartsPattern*,vector<double> > >::iterator i;
123			for (i = _vschrg.begin();i != _vschrg.end();i++)
124			if (i->first->Match(mol))
125			{
126			_mlist = i->first->GetUMapList();
127			unsigned int k;
128			vector<vector<int> >::iterator j;
129
130			for (j = _mlist.begin();j != _mlist.end();j++)
131			for (k = 0;k < j->size();k++)
132			mol.GetAtom((*j)[k])->SetPartialCharge(i->second[k]);
133			}
134			}
135
136			void OBPhModel::CorrectForPH(OBMol &mol)
137			{
138			if (!_init)
139			Init();
140			if (mol.IsCorrectedForPH())
141			return;
142			if (!mol.AutomaticFormalCharge())
143			return;
144
145			mol.SetCorrectedForPH();
146
147			obErrorLog.ThrowError(__FUNCTION__,
148			"Ran OpenBabel::CorrectForPH", obAuditMsg);
149
150			OBAtom *atom;
151			vector<OBNodeBase*>::iterator j;
152			for (atom = mol.BeginAtom(j);atom;atom = mol.NextAtom(j))
153			atom->SetFormalCharge(0);
154
155			vector<OBChemTsfm*>::iterator i;
156			for (i = _vtsfm.begin();i != _vtsfm.end();i++)
157			(*i)->Apply(mol);
158
159			atomtyper.CorrectAromaticNitrogens(mol);
160			}
161
162			// Portions of this documentation adapted from the JOELib docs, written by
163			// Joerg Wegner
164			/** \class OBChemTsfm
165			\brief SMARTS based structural modification (chemical transformation)
166
167			Transformation of chemical structures can be used for pH value correction
168			(i.e. via OBPhModel and OBMol::CorrectForPH()). The OBChemTsfm class
169			defines SMARTS based TRANSFORM patterns to delete atoms, change atom types,
170			atom formal charges, and bond types.
171
172			For storing and converting chemical reaction files, use the OBReaction class.
173			**/
174			bool OBChemTsfm::Init(string &bgn,string &end)
175			{
176			if (!_bgn.Init(bgn))
177			return(false);
178			if (!end.empty())
179			if (!_end.Init(end))
180			return(false);
181
182			//find atoms to be deleted
183			unsigned int i,j;
184			int vb;
185			bool found;
186			for (i = 0;i < _bgn.NumAtoms();i++)
187			if ((vb = _bgn.GetVectorBinding(i)))
188			{
189			found = false;
190			for (j = 0;j < _end.NumAtoms();j++)
191			if (vb == _end.GetVectorBinding(j))
192			{
193			found = true;
194			break;
195			}
196
197			if (!found)
198			_vadel.push_back(i);
199			}
200
201			//find elements to be changed
202			int ele;
203			for (i = 0;i < _bgn.NumAtoms();i++)
204			if ((vb = _bgn.GetVectorBinding(i)) != 0)
205			{
206			ele = _bgn.GetAtomicNum(i);
207			for (j = 0;j < _end.NumAtoms();j++)
208			if (vb == _end.GetVectorBinding(j))
209			if (ele != _end.GetAtomicNum(j))
210			{
211			_vele.push_back(pair<int,int> (i,_end.GetAtomicNum(j)));
212			break;
213			}
214			}
215
216			//find charges to modify
217			int chrg;
218			for (i = 0;i < _bgn.NumAtoms();i++)
219			if ((vb = _bgn.GetVectorBinding(i)))
220			{
221			chrg = _bgn.GetCharge(i);
222			for (j = 0;j < _end.NumAtoms();j++)
223			if (vb == _end.GetVectorBinding(j))
224			if (chrg != _end.GetCharge(j))
225			_vchrg.push_back(pair<int,int> (i,_end.GetCharge(j)));
226			}
227
228			//find bonds to be modified
229			//find bonds to be modified
230			int bsrc,bdst,bord,bvb1,bvb2;
231			int esrc,edst,eord,evb1,evb2;
232
233			for (i = 0;i < _bgn.NumBonds();i++)
234			{
235			_bgn.GetBond(bsrc,bdst,bord,i);
236			bvb1 = _bgn.GetVectorBinding(bsrc);
237			bvb2 = _bgn.GetVectorBinding(bdst);
238			if (!bvb1 \|\| !bvb2)
239			continue;
240
241			for (j = 0;j < _end.NumBonds();j++)
242			{
243			_end.GetBond(esrc,edst,eord,j);
244			evb1 = _end.GetVectorBinding(esrc);
245			evb2 = _end.GetVectorBinding(edst);
246			if ((bvb1 == evb1 && bvb2 == evb2) \|\| (bvb1 == evb2 && bvb2 == evb1))
247			{
248			if (bord == eord)
249			break; //nothing to modify if bond orders identical
250			_vbond.push_back(pair<pair<int,int>,int> (pair<int,int> (bsrc,bdst),eord));
251			break;
252			}
253			}
254			}
255
256			//make sure there is some kind of transform to do here
257			if (_vadel.empty() && _vchrg.empty() && _vbond.empty())
258			return(false);
259
260			return(true);
261			}
262
263			bool OBChemTsfm::Apply(OBMol &mol)
264			{
265			if (!_bgn.Match(mol))
266			return(false);
267
268			vector<vector<int> > mlist = _bgn.GetUMapList();
269
270			obErrorLog.ThrowError(__FUNCTION__,
271			"Ran OpenBabel::OBChemTransform", obAuditMsg);
272
273			if (!_vchrg.empty()) //modify charges
274			{
275			vector<vector<int> >::iterator i;
276			vector<pair<int,int> >::iterator j;
277
278			for (i = mlist.begin();i != mlist.end();i++)
279			for (j = _vchrg.begin();j != _vchrg.end();j++)
280			if (j->first < (signed)i->size()) //goof proofing
281			mol.GetAtom((*i)[j->first])->SetFormalCharge(j->second);
282
283			mol.UnsetImplicitValencePerceived();
284			}
285
286			if (!_vbond.empty()) //modify bond orders
287			{
288			OBBond *bond;
289			vector<vector<int> >::iterator i;
290			vector<pair<pair<int,int>,int> >::iterator j;
291			for (i = mlist.begin();i != mlist.end();i++)
292			for (j = _vbond.begin();j != _vbond.end();j++)
293			{
294			bond = mol.GetBond((i)[j->first.first],(i)[j->first.second]);
295			if (!bond)
296			{
297			obErrorLog.ThrowError(__FUNCTION__, "unable to find bond", obDebug);
298			continue;
299			}
300
301			bond->SetBO(j->second);
302			}
303			}
304
305			if (!_vadel.empty() \|\| !_vele.empty()) //delete atoms and change elements
306			{
307			vector<int>::iterator j;
308			vector<vector<int> >::iterator i;
309
310			if (!_vele.empty())
311			{
312			vector<pair<int,int> >::iterator k;
313			for (i = mlist.begin();i != mlist.end();i++)
314			for (k = _vele.begin();k != _vele.end();k++)
315			mol.GetAtom((*i)[k->first])->SetAtomicNum(k->second);
316			}
317
318			//make sure same atom isn't deleted twice
319			vector<bool> vda;
320			vector<OBNodeBase*> vdel;
321			vda.resize(mol.NumAtoms()+1,false);
322			for (i = mlist.begin();i != mlist.end();i++)
323			for (j = _vadel.begin();j != _vadel.end();j++)
324			if (!vda[(i)[j]])
325			{
326			vda[(i)[j]] = true;
327			vdel.push_back(mol.GetAtom((i)[j]));
328			}
329
330			vector<OBNodeBase*>::iterator k;
331			for (k = vdel.begin();k != vdel.end();k++)
332			mol.DeleteAtom((OBAtom)k);
333			}
334
335			return(true);
336			}
337
338			} //namespace OpenBabel
339
340			//! \file phmodel.cpp
341			//! \brief Read pH rules and assign charges.