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 <set>

#include "mol.hpp"
#include "obconversion.hpp"
#include "obmolecformat.hpp"

#ifdef HAVE_SSTREAM
#include <sstream>
#else
#include <strstream>
#endif

using namespace std;
namespace OpenBabel
{

class OOPSEFormat : public OBMoleculeFormat
{
public:
  //Register this format type ID
  OOPSEFormat()
  {
    OBConversion::RegisterFormat("in", this, "chemical/x-in");
  }

  virtual const char* Description() //required
  {
    return
      "oopse cartesian coordinates format\n";
  };

  virtual const char* SpecificationURL()
  { return "";}; //optional

  virtual const char* GetMIMEType() 
  { return "chemical/x-in"; };

   virtual unsigned int Flags() { return WRITEONEONLY;}

  //*** This section identical for most OBMol conversions ***
  ////////////////////////////////////////////////////
  /// The "API" interface functions
  //virtual bool ReadMolecule(OBBase* pOb, OBConversion* pConv);
  virtual bool WriteMolecule(OBBase* pOb, OBConversion* pConv);

  private:
      bool AreSameFragments(OBMol& mol, vector<int>& frag1, vector<int>& frag2);
      void findAngles(OBMol& mol);
      OBMol* createMolFromFragment(OBMol& mol, vector<int>& fragment);
      void WriteMDFile(vector<OBMol*> mols, vector<int> numMols);
      void WriteINFile(OBMol& mol, ostream& ofs, vector<int>& indices);
};


//Make an instance of the format class
OOPSEFormat theOOPSEFormat;

////////////////////////////////////////////////////////////////

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());
    }

    WriteMDFile(mdMols, numMols);    

    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, 
    for (unsigned int i =0 ; i < frag1.size(); ++i)
    {
        if (strcmp(mol.GetAtom(frag1[i])->GetType(), mol.GetAtom(frag2[i])->GetType()) )
        {
            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)
{
    std::string identLevel1("\t");
    std::string identLevel2("\t\t");
    std::string molPrefix("MolName");
    ostream& ofs = std::cout;
    const int BUFFLEN = 1024;
    char buffer[BUFFLEN];
    
    for(unsigned int i = 0; i < mols.size(); ++i)
    {
        OBMol* pmol = mols[i];
        map<OBAtom*, int> atomMap;
        ofs << "molecule {\n";
        sprintf(buffer, "%d", i);
        ofs << identLevel1 << "name = " << "\"" << molPrefix << buffer << "\"" << ";\n";

        
        //atom
        ofs << identLevel1 << "nAtoms = " << pmol->NumAtoms() << ";\n";
        int ai = 0;
        FOR_ATOMS_OF_MOL(atom, *pmol ) {
            ofs << identLevel1 << "atom[" << ai << "] {\n";
            ofs << identLevel2 << "type = " << "\"" << atom->GetType() << "\"" << ";\n";
            ofs << identLevel2 << "position(" << atom->GetX() << ", " << atom->GetY() << ", " << atom->GetZ() << ");\n";
            ofs << identLevel1 << "}\n";
            atomMap[&(*atom)] = ai++;
        }        

        //bond
        ofs << identLevel1 << "nBonds = " << pmol->NumBonds() << ";\n";
        int bi = 0;
        FOR_BONDS_OF_MOL(bond, *pmol ) {
            ofs << identLevel1 << "bond[" << bi++ << "] {\n";
            ofs << identLevel2 << "member(" << atomMap[bond->GetBeginAtom()] <<  ", " << atomMap[bond->GetEndAtom()] << ");\n";
            ofs << identLevel1 << "}\n";
        }  
        /*
        //bend
        OBGenericData* pGenericData = pmol->GetData(OBGenericDataType::AngleData);
        OBAngleData* pAngleData = dynamic_cast<OBAngleData*>(pGenericData);
        vector<OBAngle> angles = pAngleData->GetData();

        ofs << 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();
            ofs << identLevel1 << "bend[" << bendIndex++ << "] {\n";
            ofs << identLevel2 << "member(" << atomMap[bendAtoms.first] <<  ", " << atomMap[bendAtoms.second] << atomMap[bendAtoms.third] <<");\n";
            ofs << 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());            
        }

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

    ofs << "nComponents = " << mols.size() << ";\n";
    
    for(unsigned int i =0; i < mols.size(); ++i)
    {
        ofs << "component{\n";
        sprintf(buffer, "%d", i);
        ofs << "type = " << molPrefix << buffer << ";\n";
        ofs << "nMol = " << numMols[i]<< ";\n";
        ofs << "}\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,"%15s%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",
                atom->GetType(),
                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:	2440
Committed:	Wed Nov 16 19:42:11 2005 UTC (18 years, 9 months ago) by tim
File size:	10955 byte(s)
Log Message:	adding openbabel
#	User	Rev	Content
1	tim	2440	/**********************************************************************
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 <set>
17
18			#include "mol.hpp"
19			#include "obconversion.hpp"
20			#include "obmolecformat.hpp"
21
22			#ifdef HAVE_SSTREAM
23			#include <sstream>
24			#else
25			#include <strstream>
26			#endif
27
28			using namespace std;
29			namespace OpenBabel
30			{
31
32			class OOPSEFormat : public OBMoleculeFormat
33			{
34			public:
35			//Register this format type ID
36			OOPSEFormat()
37			{
38			OBConversion::RegisterFormat("in", this, "chemical/x-in");
39			}
40
41			virtual const char* Description() //required
42			{
43			return
44			"oopse cartesian coordinates format\n";
45			};
46
47			virtual const char* SpecificationURL()
48			{ return "";}; //optional
49
50			virtual const char* GetMIMEType()
51			{ return "chemical/x-in"; };
52
53			virtual unsigned int Flags() { return WRITEONEONLY;}
54
55			//* This section identical for most OBMol conversions *
56			////////////////////////////////////////////////////
57			/// The "API" interface functions
58			//virtual bool ReadMolecule(OBBase* pOb, OBConversion* pConv);
59			virtual bool WriteMolecule(OBBase* pOb, OBConversion* pConv);
60
61			private:
62			bool AreSameFragments(OBMol& mol, vector<int>& frag1, vector<int>& frag2);
63			void findAngles(OBMol& mol);
64			OBMol* createMolFromFragment(OBMol& mol, vector<int>& fragment);
65			void WriteMDFile(vector<OBMol*> mols, vector<int> numMols);
66			void WriteINFile(OBMol& mol, ostream& ofs, vector<int>& indices);
67			};
68
69
70			//Make an instance of the format class
71			OOPSEFormat theOOPSEFormat;
72
73			////////////////////////////////////////////////////////////////
74
75			bool OOPSEFormat::WriteMolecule(OBBase* pOb, OBConversion* pConv)
76			{
77			OBMol* pmol = dynamic_cast<OBMol*>(pOb);
78			if(pmol==NULL)
79			return false;
80
81			vector<vector<int> > fragmentLists;
82			pmol->ContigFragList(fragmentLists);
83			OBBitVec unused;
84			vector<bool> used(fragmentLists.size(), 0);
85			vector<vector<int> > molecules;
86			vector<int> indices;
87			for(int i =0; i < used.size(); ++i) {
88			if (used[i])
89			{
90			continue;
91			}
92			used[i] = true;
93			vector<int> sameMolTypes;
94			sameMolTypes.push_back(i);
95			indices.insert(indices.end(), fragmentLists[i].begin(), fragmentLists[i].end());
96			for (int j = i + 1;j < used.size(); ++j)
97			{
98			if (used[j])
99			{
100			continue;
101			}
102
103			if (AreSameFragments(*pmol, fragmentLists[i], fragmentLists[j]))
104			{
105			sameMolTypes.push_back(j);
106			indices.insert(indices.end(), fragmentLists[j].begin(), fragmentLists[j].end());
107			used[j]=true;
108			}
109			}
110			molecules.push_back(sameMolTypes);
111
112			}
113
114			//
115			vector<OBMol*> mdMols;
116			vector<int> numMols;
117			for(vector<vector<int> >::iterator i = molecules.begin(); i != molecules.end(); ++i)
118			{
119			mdMols.push_back(createMolFromFragment(*pmol, fragmentLists[i->front()]));
120			numMols.push_back((*i).size());
121			}
122
123			WriteMDFile(mdMols, numMols);
124
125			for(vector<OBMol*>::iterator i = mdMols.begin(); i != mdMols.end(); ++i)
126			{
127			delete *i;
128			}
129
130			//
131			WriteINFile(pmol, pConv->GetOutStream(), indices);
132
133
134
135			return(true);
136			}
137
138			bool OOPSEFormat::AreSameFragments(OBMol& mol, vector<int>& frag1, vector<int>& frag2)
139			{
140			if (frag1.size() != frag2.size())
141			{
142			return false;
143			}
144
145			//exact graph matching is a NP complete problem,
146			for (unsigned int i =0 ; i < frag1.size(); ++i)
147			{
148			if (strcmp(mol.GetAtom(frag1[i])->GetType(), mol.GetAtom(frag2[i])->GetType()) )
149			{
150			return false;
151			}
152			}
153			return true;
154			}
155
156			struct SameAngle
157			{
158			bool operator()(const triple<OBAtom,OBAtom,OBAtom> t1, const triple<OBAtom,OBAtom,OBAtom> t2) const
159			{
160			return (t1.second == t2.second) && ( (t1.first == t2.first && t1.third == t2.third) \|\| (t1.first == t2.third && t1.third == t2.first));
161			}
162			};
163
164			void OOPSEFormat::findAngles(OBMol& mol)
165			{
166			/*
167			//if already has data return
168			if(mol.HasData(OBGenericDataType::AngleData))
169			return;
170
171			vector<OBEdgeBase*>::iterator bi1,bi2;
172			OBBond* bond;
173			OBAtom a,b,*c;
174
175			set<triple<OBAtom,OBAtom,OBAtom*>, SameAngle> uniqueAngles;
176			//loop through all bonds generating torsions
177			for(bond = mol.BeginBond(bi1);bond;bond = mol.NextBond(bi1))
178			{
179			b = bond->GetBeginAtom();
180			c = bond->GetEndAtom();
181			if(b->IsHydrogen())
182			continue;
183
184			for(a = b->BeginNbrAtom(bi2);a;a = b->NextNbrAtom(bi2))
185			{
186			if(a == c)
187			continue;
188
189			uniqueAngles.insert(triple<OBAtom,OBAtom,OBAtom*>(a, b, c));
190			}
191			}
192
193			//get new data and attach it to molecule
194			OBAngleData *angles = new OBAngleData;
195			mol.SetData(angles);
196			set<triple<OBAtom,OBAtom,OBAtom*>, SameAngle>::iterator i;
197
198			for (i = uniqueAngles.begin(); i != uniqueAngles.end(); ++i) {
199			OBAngle angle;
200			angle.SetAtoms(i->first, i->second, i->second);
201			angles->SetData(angle);
202			}
203			*/
204			}
205			OBMol* OOPSEFormat::createMolFromFragment(OBMol& mol, vector<int>& fragment)
206			{
207			OBMol* newMol = new OBMol();
208			newMol->ReserveAtoms(fragment.size());
209			newMol->BeginModify();
210			for(vector<int>::iterator i = fragment.begin(); i != fragment.end(); ++i)
211			{
212			OBAtom* newAtom = newMol->NewAtom();
213			newAtom = mol.GetAtom(*i);
214			}
215			newMol->EndModify();
216			newMol->ConnectTheDots();
217			findAngles(*newMol);
218			newMol->FindTorsions();
219			return newMol;
220			}
221			void OOPSEFormat::WriteMDFile(vector<OBMol*> mols, vector<int> numMols)
222			{
223			std::string identLevel1("\t");
224			std::string identLevel2("\t\t");
225			std::string molPrefix("MolName");
226			ostream& ofs = std::cout;
227			const int BUFFLEN = 1024;
228			char buffer[BUFFLEN];
229
230			for(unsigned int i = 0; i < mols.size(); ++i)
231			{
232			OBMol* pmol = mols[i];
233			map<OBAtom*, int> atomMap;
234			ofs << "molecule {\n";
235			sprintf(buffer, "%d", i);
236			ofs << identLevel1 << "name = " << "\"" << molPrefix << buffer << "\"" << ";\n";
237
238
239			//atom
240			ofs << identLevel1 << "nAtoms = " << pmol->NumAtoms() << ";\n";
241			int ai = 0;
242			FOR_ATOMS_OF_MOL(atom, *pmol ) {
243			ofs << identLevel1 << "atom[" << ai << "] {\n";
244			ofs << identLevel2 << "type = " << "\"" << atom->GetType() << "\"" << ";\n";
245			ofs << identLevel2 << "position(" << atom->GetX() << ", " << atom->GetY() << ", " << atom->GetZ() << ");\n";
246			ofs << identLevel1 << "}\n";
247			atomMap[&(*atom)] = ai++;
248			}
249
250			//bond
251			ofs << identLevel1 << "nBonds = " << pmol->NumBonds() << ";\n";
252			int bi = 0;
253			FOR_BONDS_OF_MOL(bond, *pmol ) {
254			ofs << identLevel1 << "bond[" << bi++ << "] {\n";
255			ofs << identLevel2 << "member(" << atomMap[bond->GetBeginAtom()] << ", " << atomMap[bond->GetEndAtom()] << ");\n";
256			ofs << identLevel1 << "}\n";
257			}
258			/*
259			//bend
260			OBGenericData* pGenericData = pmol->GetData(OBGenericDataType::AngleData);
261			OBAngleData* pAngleData = dynamic_cast<OBAngleData*>(pGenericData);
262			vector<OBAngle> angles = pAngleData->GetData();
263
264			ofs << identLevel1 << "nBends = " << angles.size() << ";\n";
265			int bendIndex = 0;
266			for (vector<OBAngle>::iterator ti = angles.begin(); ti != angles.end(); ++ti)
267			{
268			triple<OBAtom, OBAtom, OBAtom*> bendAtoms = ti->getAtoms();
269			ofs << identLevel1 << "bend[" << bendIndex++ << "] {\n";
270			ofs << identLevel2 << "member(" << atomMap[bendAtoms.first] << ", " << atomMap[bendAtoms.second] << atomMap[bendAtoms.third] <<");\n";
271			ofs << identLevel1 << "}\n";
272			}
273
274			//torsion
275			pGenericData = pmol->GetData(OBGenericDataType::TorsionData);
276			OBTorsionData* pTorsionData = dynamic_cast<OBTorsionData*>(pGenericData);
277			vector<OBTorsion> torsions = pTorsionData->GetData();
278			vector<quad<OBAtom,OBAtom,OBAtom,OBAtom> > torsionArray;
279			for (vector<OBTorsion>::iterator ti = torsions.begin(); ti != torsions.end(); ++ti)
280			{
281			vector<quad<OBAtom,OBAtom,OBAtom,OBAtom> > tmpTorsions = ti->getTorsions();
282			torsionArray.insert(torsionArray.end(), tmpTorsions.begin(), tmpTorsions.end());
283			}
284
285			ofs << identLevel1 << "nTorsions = " << torsionArray.size() << ";\n";
286			int torsionIndex = 0;
287			for (vector<quad<OBAtom,OBAtom,OBAtom,OBAtom> >::iterator ti = torsionArray.begin(); ti != torsionArray.end(); ++ti)
288			{
289			ofs << identLevel1 << "torsion[" << torsionIndex++ << "] {\n";
290			ofs << identLevel2 << "member(" << atomMap[ti->first] << ", " << atomMap[ti->second] <<", " << atomMap[ti->third] <<", " << atomMap[ti->forth] << ");\n";
291			ofs << identLevel1 << "}\n";
292			}
293			*/
294			ofs << "}\n";
295			}
296
297			ofs << "nComponents = " << mols.size() << ";\n";
298
299			for(unsigned int i =0; i < mols.size(); ++i)
300			{
301			ofs << "component{\n";
302			sprintf(buffer, "%d", i);
303			ofs << "type = " << molPrefix << buffer << ";\n";
304			ofs << "nMol = " << numMols[i]<< ";\n";
305			ofs << "}\n";
306			}
307			}
308			void OOPSEFormat::WriteINFile(OBMol& mol, ostream& ofs, vector<int>& indices)
309			{
310			unsigned int i;
311			char buffer[BUFF_SIZE];
312
313			sprintf(buffer,"%d", mol.NumAtoms());
314			ofs << buffer << endl;
315			//sprintf(buffer,"0;%f, 0, 0; 0, %15.7f, 0; 0, 0, %15.7f;", boxx, boxy, boxz);
316			sprintf(buffer, "0;%f, 0, 0; 0, %15.7f, 0; 0, 0, %15.7f;", 100.0, 100.0, 100.0);
317			ofs << buffer << endl;
318
319			OBAtom *atom;
320			string str,str1;
321
322			for(vector<int>::iterator i = indices.begin();i != indices.end(); ++i)
323			{
324			atom = mol.GetAtom(*i);
325			sprintf(buffer,"%15s%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",
326			atom->GetType(),
327			atom->GetX(), atom->GetY(), atom->GetZ(),
328			0.0, 0.0, 0.0,
329			0.0, 0.0, 0.0, 0.0,
330			0.0, 0.0, 0.0
331			);
332			ofs << buffer << endl;
333			}
334
335			}
336
337			} //namespace OpenBabel
338