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.
#	Content
1	/**********************************************************************
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	STR_DEFINE(_dir, FRC_PATH);
41	_envvar = "FORCE_PARAM_PATH";
42	_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.