src/openbabel/oopseformat.cpp

/**********************************************************************
Copyright (C) 2000 by OpenEye Scientific Software, Inc.
Some portions Copyright (C) 2001-2005 by Geoffrey R. Hutchison
Some portions Copyright (C) 2004 by Chris Morley
 
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 "oopseformat.hpp"
#include <fstream>
namespace OpenBabel
{

bool OOPSEFormat::WriteMolecule(OBBase* pOb, OBConversion* pConv)
{
    OBMol* pmol = dynamic_cast<OBMol*>(pOb);
    if(pmol==NULL)
        return false;
    
    vector<vector<int> > fragmentLists;
    pmol->ContigFragList(fragmentLists);
    OBBitVec unused;
    vector<bool> used(fragmentLists.size(), 0);
    vector<vector<int> > molecules;
    vector<int> indices;
    for(int i =0; i < used.size(); ++i) {
        if (used[i])
        {
            continue;
        }
        used[i] = true;
        vector<int> sameMolTypes;
        sameMolTypes.push_back(i);
        indices.insert(indices.end(), fragmentLists[i].begin(), fragmentLists[i].end());
        for (int j = i + 1;j < used.size(); ++j)
        {
            if (used[j])
            {
                continue;
            }
            
            if (AreSameFragments(*pmol, fragmentLists[i], fragmentLists[j]))
            {
                sameMolTypes.push_back(j);
                indices.insert(indices.end(), fragmentLists[j].begin(), fragmentLists[j].end());
                used[j]=true;
            }
        }
        molecules.push_back(sameMolTypes);
        
    }

    //
    vector<OBMol*> mdMols;    
    vector<int> numMols;
    for(vector<vector<int> >::iterator  i = molecules.begin(); i != molecules.end(); ++i) 
    {
        mdMols.push_back(createMolFromFragment(*pmol, fragmentLists[i->front()]));
        numMols.push_back((*i).size());
    }

    string OutputFileName = pConv->GetInFilename();
    unsigned int pos = OutputFileName.rfind(".");
    if(pos==string::npos)
        OutputFileName += ".md";
    else
        OutputFileName = OutputFileName.substr(0, pos) + ".md";       
    ofstream ofs(OutputFileName.c_str());
    if(!ofs)
    {
         cerr << "Cannot write to " << OutputFileName <<endl;
         return false;
    }
                
    WriteMDFile(mdMols, numMols, ofs);    

    for(vector<OBMol*>::iterator  i = mdMols.begin(); i != mdMols.end(); ++i) 
    {
        delete *i;
    }

    //    
    WriteINFile(*pmol, *pConv->GetOutStream(), indices);


    return(true);
}

bool OOPSEFormat::AreSameFragments(OBMol& mol, vector<int>& frag1, vector<int>& frag2)
{
    if (frag1.size() != frag2.size())
    {
        return false;
    }

    //exact graph matching is a NP complete problem
    /** @todo using sparse matrix to store the connectivities*/
    for (unsigned int i =0 ; i < frag1.size(); ++i)
    {
        OBAtom* atom1 = mol.GetAtom(frag1[i]);
        OBAtom* atom2 = mol.GetAtom(frag2[i]);
        
        if (atom1->GetAtomicNum() != atom2->GetAtomicNum())
        {
            return false;
        }
    }
    return true;
}

struct SameAngle
{
  bool operator()(const triple<OBAtom*,OBAtom*,OBAtom*> t1, const triple<OBAtom*,OBAtom*,OBAtom*> t2) const
  {
    return (t1.second == t2.second) && ( (t1.first == t2.first && t1.third == t2.third) || (t1.first == t2.third && t1.third == t2.first));
  }
};

void OOPSEFormat::findAngles(OBMol& mol)
{
    /*
    //if already has data return
    if(mol.HasData(OBGenericDataType::AngleData))
        return;

    vector<OBEdgeBase*>::iterator bi1,bi2;
    OBBond* bond;
    OBAtom *a,*b,*c;

    set<triple<OBAtom*,OBAtom*,OBAtom*>, SameAngle> uniqueAngles;
    //loop through all bonds generating torsions
    for(bond = mol.BeginBond(bi1);bond;bond = mol.NextBond(bi1))
    {
        b = bond->GetBeginAtom();
        c = bond->GetEndAtom();
        if(b->IsHydrogen())
            continue;

        for(a = b->BeginNbrAtom(bi2);a;a = b->NextNbrAtom(bi2))
        {
            if(a == c)
                continue;          
            
            uniqueAngles.insert(triple<OBAtom*,OBAtom*,OBAtom*>(a, b, c));
        }
    }

    //get new data and attach it to molecule
    OBAngleData *angles = new OBAngleData;
    mol.SetData(angles);
    set<triple<OBAtom*,OBAtom*,OBAtom*>, SameAngle>::iterator i;

    for (i = uniqueAngles.begin(); i != uniqueAngles.end(); ++i) {
        OBAngle angle;
        angle.SetAtoms(i->first, i->second, i->second);
        angles->SetData(angle);
    }
    */
}
OBMol* OOPSEFormat::createMolFromFragment(OBMol& mol, vector<int>& fragment)
{
    OBMol* newMol = new OBMol();
    newMol->ReserveAtoms(fragment.size());
    newMol->BeginModify();
    for(vector<int>::iterator i = fragment.begin(); i != fragment.end(); ++i)
    {
        OBAtom* newAtom = newMol->NewAtom();
        *newAtom = *mol.GetAtom(*i);
    }
    newMol->EndModify();
    newMol->ConnectTheDots();
    findAngles(*newMol);
    newMol->FindTorsions();
    return newMol;
}
void OOPSEFormat::WriteMDFile(vector<OBMol*> mols, vector<int> numMols, ostream& os)
{
    std::string identLevel1("\t");
    std::string identLevel2("\t\t");
    std::string molPrefix("MolName");
    const int BUFFLEN = 1024;
    char buffer[BUFFLEN];
    
    for(unsigned int i = 0; i < mols.size(); ++i)
    {
        OBMol* pmol = mols[i];
        map<OBAtom*, int> atomMap;
        os << "molecule {\n";
        sprintf(buffer, "%d", i);
        os << identLevel1 << "name = " << "\"" << molPrefix << buffer << "\"" << ";\n";

        
        //atom
        int ai = 0;
        FOR_ATOMS_OF_MOL(atom, *pmol ) {
            os << identLevel1 << "atom[" << ai << "] {\n";
            os << identLevel2 << "type = " << "\"" << etab.GetSymbol(atom->GetAtomicNum()) << "\"" << ";\n";
            os << identLevel1 << "}\n";
            atomMap[&(*atom)] = ai++;
        }        
        os << "\n";

        //bond
        FOR_BONDS_OF_MOL(bond, *pmol ) {
            os << identLevel1 << "bond {\n";
            os << identLevel2 << "members(" << atomMap[bond->GetBeginAtom()] <<  ", " << atomMap[bond->GetEndAtom()] << ");\n";
            os << identLevel1 << "}\n";
        }  
        /*
        //bend
        OBGenericData* pGenericData = pmol->GetData(OBGenericDataType::AngleData);
        OBAngleData* pAngleData = dynamic_cast<OBAngleData*>(pGenericData);
        vector<OBAngle> angles = pAngleData->GetData();

        os << identLevel1 << "nBends = " << angles.size() << ";\n";        
        int bendIndex = 0;
        for (vector<OBAngle>::iterator ti = angles.begin(); ti != angles.end(); ++ti)
        {
            triple<OBAtom*, OBAtom*, OBAtom*> bendAtoms = ti->getAtoms();
            os << identLevel1 << "bend[" << bendIndex++ << "] {\n";
            os << identLevel2 << "member(" << atomMap[bendAtoms.first] <<  ", " << atomMap[bendAtoms.second] << atomMap[bendAtoms.third] <<");\n";
            os << identLevel1 << "}\n";            
        }
        
        //torsion
        pGenericData = pmol->GetData(OBGenericDataType::TorsionData);
        OBTorsionData* pTorsionData = dynamic_cast<OBTorsionData*>(pGenericData);
        vector<OBTorsion> torsions = pTorsionData->GetData();
        vector<quad<OBAtom*,OBAtom*,OBAtom*,OBAtom*> > torsionArray;
        for (vector<OBTorsion>::iterator ti = torsions.begin(); ti != torsions.end(); ++ti)
        {
            vector<quad<OBAtom*,OBAtom*,OBAtom*,OBAtom*> > tmpTorsions = ti->getTorsions();
            torsionArray.insert(torsionArray.end(), tmpTorsions.begin(), tmpTorsions.end());            
        }

        os << identLevel1 << "nTorsions = " << torsionArray.size() << ";\n";
        int torsionIndex = 0;
        for (vector<quad<OBAtom*,OBAtom*,OBAtom*,OBAtom*> >::iterator ti = torsionArray.begin(); ti != torsionArray.end(); ++ti)
        {
            os << identLevel1 << "torsion[" << torsionIndex++ << "] {\n";
            os << identLevel2 << "member(" << atomMap[ti->first] <<  ", " << atomMap[ti->second] <<", " << atomMap[ti->third] <<", " << atomMap[ti->forth] << ");\n";
            os << identLevel1 << "}\n";          
        }
        */
        os << "}\n";
        os << "\n";

    }

    os << "\n";
    os << "nComponents = " << mols.size() << ";\n";
    
    for(unsigned int i =0; i < mols.size(); ++i)
    {
        os << "component{\n";
        sprintf(buffer, "%d", i);
        os << "type = " << molPrefix << buffer << ";\n";
        os << "nMol = " << numMols[i]<< ";\n";
        os << "}\n";
    }
}
void OOPSEFormat::WriteINFile(OBMol& mol, ostream& ofs, vector<int>& indices)
{
    unsigned int i;
    char buffer[BUFF_SIZE];
    
    sprintf(buffer,"%d", mol.NumAtoms());
    ofs << buffer << endl;
    //sprintf(buffer,"0;%f, 0, 0; 0, %15.7f, 0; 0, 0, %15.7f;", boxx, boxy, boxz);
    sprintf(buffer, "0;%f, 0, 0; 0, %15.7f, 0; 0, 0, %15.7f;", 100.0, 100.0, 100.0);
    ofs << buffer << endl;

    OBAtom *atom;
    string str,str1;
    
    for(vector<int>::iterator i = indices.begin();i != indices.end(); ++i)
    {
        atom = mol.GetAtom(*i);
        sprintf(buffer,"%-10s%-15.5f%-15.5f%-15.5f%-15.5f%-15.5f%-15.5f%-15.5f%-15.5f%-15.5f%-15.5f%-15.5f%-15.5f%-15.5f",
                etab.GetSymbol(atom->GetAtomicNum()),
                atom->GetX(), atom->GetY(), atom->GetZ(),
                0.0, 0.0, 0.0,
                0.0, 0.0, 0.0, 0.0,
                0.0, 0.0, 0.0
                );
        ofs << buffer << endl;
    }
    
}

} //namespace OpenBabel

Revision:	2469
Committed:	Fri Dec 2 15:38:03 2005 UTC (18 years, 7 months ago) by tim
File size:	9888 byte(s)
Log Message:	End of the Link --> List Return of the Oject-Oriented replace yacc/lex parser with antlr parser
#	Content
1	/**********************************************************************
2	Copyright (C) 2000 by OpenEye Scientific Software, Inc.
3	Some portions Copyright (C) 2001-2005 by Geoffrey R. Hutchison
4	Some portions Copyright (C) 2004 by Chris Morley
5
6	This program is free software; you can redistribute it and/or modify
7	it under the terms of the GNU General Public License as published by
8	the Free Software Foundation version 2 of the License.
9
10	This program is distributed in the hope that it will be useful,
11	but WITHOUT ANY WARRANTY; without even the implied warranty of
12	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13	GNU General Public License for more details.
14	***********************************************************************/
15
16	#include "oopseformat.hpp"
17	#include <fstream>
18	namespace OpenBabel
19	{
20
21	bool OOPSEFormat::WriteMolecule(OBBase* pOb, OBConversion* pConv)
22	{
23	OBMol* pmol = dynamic_cast<OBMol*>(pOb);
24	if(pmol==NULL)
25	return false;
26
27	vector<vector<int> > fragmentLists;
28	pmol->ContigFragList(fragmentLists);
29	OBBitVec unused;
30	vector<bool> used(fragmentLists.size(), 0);
31	vector<vector<int> > molecules;
32	vector<int> indices;
33	for(int i =0; i < used.size(); ++i) {
34	if (used[i])
35	{
36	continue;
37	}
38	used[i] = true;
39	vector<int> sameMolTypes;
40	sameMolTypes.push_back(i);
41	indices.insert(indices.end(), fragmentLists[i].begin(), fragmentLists[i].end());
42	for (int j = i + 1;j < used.size(); ++j)
43	{
44	if (used[j])
45	{
46	continue;
47	}
48
49	if (AreSameFragments(*pmol, fragmentLists[i], fragmentLists[j]))
50	{
51	sameMolTypes.push_back(j);
52	indices.insert(indices.end(), fragmentLists[j].begin(), fragmentLists[j].end());
53	used[j]=true;
54	}
55	}
56	molecules.push_back(sameMolTypes);
57
58	}
59
60	//
61	vector<OBMol*> mdMols;
62	vector<int> numMols;
63	for(vector<vector<int> >::iterator i = molecules.begin(); i != molecules.end(); ++i)
64	{
65	mdMols.push_back(createMolFromFragment(*pmol, fragmentLists[i->front()]));
66	numMols.push_back((*i).size());
67	}
68
69	string OutputFileName = pConv->GetInFilename();
70	unsigned int pos = OutputFileName.rfind(".");
71	if(pos==string::npos)
72	OutputFileName += ".md";
73	else
74	OutputFileName = OutputFileName.substr(0, pos) + ".md";
75	ofstream ofs(OutputFileName.c_str());
76	if(!ofs)
77	{
78	cerr << "Cannot write to " << OutputFileName <<endl;
79	return false;
80	}
81
82	WriteMDFile(mdMols, numMols, ofs);
83
84	for(vector<OBMol*>::iterator i = mdMols.begin(); i != mdMols.end(); ++i)
85	{
86	delete *i;
87	}
88
89	//
90	WriteINFile(pmol, pConv->GetOutStream(), indices);
91
92
93
94	return(true);
95	}
96
97	bool OOPSEFormat::AreSameFragments(OBMol& mol, vector<int>& frag1, vector<int>& frag2)
98	{
99	if (frag1.size() != frag2.size())
100	{
101	return false;
102	}
103
104	//exact graph matching is a NP complete problem
105	/** @todo using sparse matrix to store the connectivities*/
106	for (unsigned int i =0 ; i < frag1.size(); ++i)
107	{
108	OBAtom* atom1 = mol.GetAtom(frag1[i]);
109	OBAtom* atom2 = mol.GetAtom(frag2[i]);
110
111	if (atom1->GetAtomicNum() != atom2->GetAtomicNum())
112	{
113	return false;
114	}
115	}
116	return true;
117	}
118
119	struct SameAngle
120	{
121	bool operator()(const triple<OBAtom,OBAtom,OBAtom> t1, const triple<OBAtom,OBAtom,OBAtom> t2) const
122	{
123	return (t1.second == t2.second) && ( (t1.first == t2.first && t1.third == t2.third) \|\| (t1.first == t2.third && t1.third == t2.first));
124	}
125	};
126
127	void OOPSEFormat::findAngles(OBMol& mol)
128	{
129	/*
130	//if already has data return
131	if(mol.HasData(OBGenericDataType::AngleData))
132	return;
133
134	vector<OBEdgeBase*>::iterator bi1,bi2;
135	OBBond* bond;
136	OBAtom a,b,*c;
137
138	set<triple<OBAtom,OBAtom,OBAtom*>, SameAngle> uniqueAngles;
139	//loop through all bonds generating torsions
140	for(bond = mol.BeginBond(bi1);bond;bond = mol.NextBond(bi1))
141	{
142	b = bond->GetBeginAtom();
143	c = bond->GetEndAtom();
144	if(b->IsHydrogen())
145	continue;
146
147	for(a = b->BeginNbrAtom(bi2);a;a = b->NextNbrAtom(bi2))
148	{
149	if(a == c)
150	continue;
151
152	uniqueAngles.insert(triple<OBAtom,OBAtom,OBAtom*>(a, b, c));
153	}
154	}
155
156	//get new data and attach it to molecule
157	OBAngleData *angles = new OBAngleData;
158	mol.SetData(angles);
159	set<triple<OBAtom,OBAtom,OBAtom*>, SameAngle>::iterator i;
160
161	for (i = uniqueAngles.begin(); i != uniqueAngles.end(); ++i) {
162	OBAngle angle;
163	angle.SetAtoms(i->first, i->second, i->second);
164	angles->SetData(angle);
165	}
166	*/
167	}
168	OBMol* OOPSEFormat::createMolFromFragment(OBMol& mol, vector<int>& fragment)
169	{
170	OBMol* newMol = new OBMol();
171	newMol->ReserveAtoms(fragment.size());
172	newMol->BeginModify();
173	for(vector<int>::iterator i = fragment.begin(); i != fragment.end(); ++i)
174	{
175	OBAtom* newAtom = newMol->NewAtom();
176	newAtom = mol.GetAtom(*i);
177	}
178	newMol->EndModify();
179	newMol->ConnectTheDots();
180	findAngles(*newMol);
181	newMol->FindTorsions();
182	return newMol;
183	}
184	void OOPSEFormat::WriteMDFile(vector<OBMol*> mols, vector<int> numMols, ostream& os)
185	{
186	std::string identLevel1("\t");
187	std::string identLevel2("\t\t");
188	std::string molPrefix("MolName");
189	const int BUFFLEN = 1024;
190	char buffer[BUFFLEN];
191
192	for(unsigned int i = 0; i < mols.size(); ++i)
193	{
194	OBMol* pmol = mols[i];
195	map<OBAtom*, int> atomMap;
196	os << "molecule {\n";
197	sprintf(buffer, "%d", i);
198	os << identLevel1 << "name = " << "\"" << molPrefix << buffer << "\"" << ";\n";
199
200
201	//atom
202	int ai = 0;
203	FOR_ATOMS_OF_MOL(atom, *pmol ) {
204	os << identLevel1 << "atom[" << ai << "] {\n";
205	os << identLevel2 << "type = " << "\"" << etab.GetSymbol(atom->GetAtomicNum()) << "\"" << ";\n";
206	os << identLevel1 << "}\n";
207	atomMap[&(*atom)] = ai++;
208	}
209	os << "\n";
210
211	//bond
212	FOR_BONDS_OF_MOL(bond, *pmol ) {
213	os << identLevel1 << "bond {\n";
214	os << identLevel2 << "members(" << atomMap[bond->GetBeginAtom()] << ", " << atomMap[bond->GetEndAtom()] << ");\n";
215	os << identLevel1 << "}\n";
216	}
217	/*
218	//bend
219	OBGenericData* pGenericData = pmol->GetData(OBGenericDataType::AngleData);
220	OBAngleData* pAngleData = dynamic_cast<OBAngleData*>(pGenericData);
221	vector<OBAngle> angles = pAngleData->GetData();
222
223	os << identLevel1 << "nBends = " << angles.size() << ";\n";
224	int bendIndex = 0;
225	for (vector<OBAngle>::iterator ti = angles.begin(); ti != angles.end(); ++ti)
226	{
227	triple<OBAtom, OBAtom, OBAtom*> bendAtoms = ti->getAtoms();
228	os << identLevel1 << "bend[" << bendIndex++ << "] {\n";
229	os << identLevel2 << "member(" << atomMap[bendAtoms.first] << ", " << atomMap[bendAtoms.second] << atomMap[bendAtoms.third] <<");\n";
230	os << identLevel1 << "}\n";
231	}
232
233	//torsion
234	pGenericData = pmol->GetData(OBGenericDataType::TorsionData);
235	OBTorsionData* pTorsionData = dynamic_cast<OBTorsionData*>(pGenericData);
236	vector<OBTorsion> torsions = pTorsionData->GetData();
237	vector<quad<OBAtom,OBAtom,OBAtom,OBAtom> > torsionArray;
238	for (vector<OBTorsion>::iterator ti = torsions.begin(); ti != torsions.end(); ++ti)
239	{
240	vector<quad<OBAtom,OBAtom,OBAtom,OBAtom> > tmpTorsions = ti->getTorsions();
241	torsionArray.insert(torsionArray.end(), tmpTorsions.begin(), tmpTorsions.end());
242	}
243
244	os << identLevel1 << "nTorsions = " << torsionArray.size() << ";\n";
245	int torsionIndex = 0;
246	for (vector<quad<OBAtom,OBAtom,OBAtom,OBAtom> >::iterator ti = torsionArray.begin(); ti != torsionArray.end(); ++ti)
247	{
248	os << identLevel1 << "torsion[" << torsionIndex++ << "] {\n";
249	os << identLevel2 << "member(" << atomMap[ti->first] << ", " << atomMap[ti->second] <<", " << atomMap[ti->third] <<", " << atomMap[ti->forth] << ");\n";
250	os << identLevel1 << "}\n";
251	}
252	*/
253	os << "}\n";
254	os << "\n";
255
256	}
257
258	os << "\n";
259	os << "nComponents = " << mols.size() << ";\n";
260
261	for(unsigned int i =0; i < mols.size(); ++i)
262	{
263	os << "component{\n";
264	sprintf(buffer, "%d", i);
265	os << "type = " << molPrefix << buffer << ";\n";
266	os << "nMol = " << numMols[i]<< ";\n";
267	os << "}\n";
268	}
269	}
270	void OOPSEFormat::WriteINFile(OBMol& mol, ostream& ofs, vector<int>& indices)
271	{
272	unsigned int i;
273	char buffer[BUFF_SIZE];
274
275	sprintf(buffer,"%d", mol.NumAtoms());
276	ofs << buffer << endl;
277	//sprintf(buffer,"0;%f, 0, 0; 0, %15.7f, 0; 0, 0, %15.7f;", boxx, boxy, boxz);
278	sprintf(buffer, "0;%f, 0, 0; 0, %15.7f, 0; 0, 0, %15.7f;", 100.0, 100.0, 100.0);
279	ofs << buffer << endl;
280
281	OBAtom *atom;
282	string str,str1;
283
284	for(vector<int>::iterator i = indices.begin();i != indices.end(); ++i)
285	{
286	atom = mol.GetAtom(*i);
287	sprintf(buffer,"%-10s%-15.5f%-15.5f%-15.5f%-15.5f%-15.5f%-15.5f%-15.5f%-15.5f%-15.5f%-15.5f%-15.5f%-15.5f%-15.5f",
288	etab.GetSymbol(atom->GetAtomicNum()),
289	atom->GetX(), atom->GetY(), atom->GetZ(),
290	0.0, 0.0, 0.0,
291	0.0, 0.0, 0.0, 0.0,
292	0.0, 0.0, 0.0
293	);
294	ofs << buffer << endl;
295	}
296
297	}
298
299	} //namespace OpenBabel
300