applications/atom2md/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
Some portions Copyright (C) 2008 by J. Daniel Gezelter
 
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 <openbabel/babelconfig.h>
#include <openbabel/obmolecformat.h>
#include <fstream>

using namespace std;
namespace OpenBabel
{
  
  class OOPSEFormat : public OBMoleculeFormat
  {
  public:
    //Register this format type ID
    OOPSEFormat() 
    {      
      OBConversion::RegisterFormat("md",this);
    }
    
    virtual const char* Description() //required
    {
      return
        "OOPSE combined meta-data / cartesian coordinates format\nNo comments yet\n";
    };
    
    virtual const char* SpecificationURL()
    {return "http://www.oopse.org";}; //optional
    
    virtual const char* GetMIMEType() 
    {return "chemical/x-md"; };
    
    virtual unsigned int Flags() 
    { 
      return NOTREADABLE | WRITEONEONLY;
    }
    
    //*** This section identical for most OBMol conversions ***
    ////////////////////////////////////////////////////
    /// The "API" interface functions
    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, ostream& os, OBMol& mol, 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());
      }
    
    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, *pmol, indices);
    
    for(vector<OBMol*>::iterator  i = mdMols.begin(); i != mdMols.end(); ++i) 
      {
        delete *i;
      }
    
    //    
    
    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, OBMol& mol, vector<int>& indices) {
    std::string indentLevel1("  ");
    std::string indentLevel2("    ");
    std::string molPrefix("MolName");
    unsigned int i;
    const int BUFFLEN = 1024;
    char buffer[BUFFLEN];
    string str, str1;
    
    
    os << "<OOPSE version=4>" << endl;
    os << "  <MetaData>" << endl << endl;
    
    for(i = 0; i < mols.size(); ++i) {
      OBMol* pmol = mols[i];
      
      pmol->ConnectTheDots();
      pmol->PerceiveBondOrders();
      //pmol->FindSSSR();
      //pmol->SetAromaticPerceived();
      //pmol->Kekulize();
      
      map<OBAtom*, int> atomMap;
      os << "molecule {\n";
      sprintf(buffer, "%d", i);
      os << indentLevel1 << "name = " << "\"" << molPrefix << buffer << "\"" << ";\n";
      
      //atom
      int ai = 0;
      FOR_ATOMS_OF_MOL(atom, *pmol ) {
        str = atom->GetType();
        ttab.SetFromType("INT");
        ttab.SetToType("INT");
        ttab.Translate(str1,str);
        os << indentLevel1 << "atom[" << ai << "] {\n";
        // os << indentLevel2 << "type = " << "\"" << etab.GetSymbol(atom->GetAtomicNum()) << "\"" << ";\n";
        os << indentLevel2 << "type = " << "\"" << str1 << "\"" << ";\n";
        os << indentLevel1 << "}\n";
        atomMap[&(*atom)] = ai++;
      }        
      os << "\n";
      
      //bond
      FOR_BONDS_OF_MOL(bond, *pmol ) {
        os << indentLevel1 << "bond {\n";
        os << indentLevel2 << "members(" << atomMap[bond->GetBeginAtom()] <<  ", " << atomMap[bond->GetEndAtom()] << ");\n";
        os << indentLevel1 << "}\n";
      }  
      /*
      //bend
      OBGenericData* pGenericData = pmol->GetData(OBGenericDataType::AngleData);
      OBAngleData* pAngleData = dynamic_cast<OBAngleData*>(pGenericData);
      vector<OBAngle> angles = pAngleData->GetData();
      
      os << indentLevel1 << "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 << indentLevel1 << "bend[" << bendIndex++ << "] {\n";
          os << indentLevel2 << "member(" << atomMap[bendAtoms.first] <<  ", " << atomMap[bendAtoms.second] << atomMap[bendAtoms.third] <<");\n";
          os << indentLevel1 << "}\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 << indentLevel1 << "nTorsions = " << torsionArray.size() << ";\n";
          int torsionIndex = 0;
          for (vector<quad<OBAtom*,OBAtom*,OBAtom*,OBAtom*> >::iterator ti = torsionArray.begin(); ti != torsionArray.end(); ++ti)
          {
          os << indentLevel1 << "torsion[" << torsionIndex++ << "] {\n";
          os << indentLevel2 << "member(" << atomMap[ti->first] <<  ", " << atomMap[ti->second] <<", " << atomMap[ti->third] <<", " << atomMap[ti->forth] << ");\n";
          os << indentLevel1 << "}\n";          
          }
      */
      os << "}" << endl;
      os << endl;
      
    }
    
    os << endl;
    
    
    for(i=0; i < mols.size(); ++i) {      
      os << "component{" << endl;
      sprintf(buffer, "%d", i);
      os << indentLevel1 << "type = " << molPrefix << buffer << ";" << endl;
      os << indentLevel1 << "nMol = " << numMols[i]<< ";" << endl;
      os << "}" << endl;
    }
    
    os << "  </MetaData>" << endl;
    os << "  <Snapshot>" << endl;
    os << "    <FrameData>" << endl;
    
    sprintf(buffer, "        Time: %.10g", 0.0);
    
    os << buffer << endl;
    
    sprintf(buffer, "        Hmat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}", 100.0, 0.0, 0.0, 0.0, 100.0, 0.0, 0.0, 0.0, 100.0);
    
    os << buffer << endl;
    os << "    </FrameData>" << endl;
    os << "    <StuntDoubles>" << endl;
    
    OBAtom *atom;
    
    for(vector<int>::iterator i = indices.begin();i != indices.end(); ++i) {     
      atom = mol.GetAtom(*i);
      sprintf(buffer, "%10d %7s %18.10g %18.10g %18.10g %13e %13e %13e", *i - 1, "pv", atom->GetX(), atom->GetY(), atom->GetZ(), 0.0, 0.0, 0.0);
      os << buffer << endl;
    }
    os << "    </StuntDoubles>" << endl;
    os << "  </Snapshot>" << endl;
    os << "</OOPSE>" << endl;
  }
  
} //namespace OpenBabel

Revision:	1210
Committed:	Wed Jan 23 03:45:33 2008 UTC (17 years, 3 months ago) by gezelter
File size:	11923 byte(s)
Log Message:	Removed older version of openbabel from our code. We now have a configure check to see if openbabel is installed and then we link to the stuff we need. Conversion to OOPSE's md format is handled by only one application (atom2md), so most of the work went on there. ElementsTable still needs some work to function in parallel.
#	User	Rev	Content
1	gezelter	1210	/**********************************************************************
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			Some portions Copyright (C) 2008 by J. Daniel Gezelter
6
7			This program is free software; you can redistribute it and/or modify
8			it under the terms of the GNU General Public License as published by
9			the Free Software Foundation version 2 of the License.
10
11			This program is distributed in the hope that it will be useful,
12			but WITHOUT ANY WARRANTY; without even the implied warranty of
13			MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14			GNU General Public License for more details.
15			***********************************************************************/
16
17			#include <openbabel/babelconfig.h>
18			#include <openbabel/obmolecformat.h>
19			#include <fstream>
20
21			using namespace std;
22			namespace OpenBabel
23			{
24
25			class OOPSEFormat : public OBMoleculeFormat
26			{
27			public:
28			//Register this format type ID
29			OOPSEFormat()
30			{
31			OBConversion::RegisterFormat("md",this);
32			}
33
34			virtual const char* Description() //required
35			{
36			return
37			"OOPSE combined meta-data / cartesian coordinates format\nNo comments yet\n";
38			};
39
40			virtual const char* SpecificationURL()
41			{return "http://www.oopse.org";}; //optional
42
43			virtual const char* GetMIMEType()
44			{return "chemical/x-md"; };
45
46			virtual unsigned int Flags()
47			{
48			return NOTREADABLE \| WRITEONEONLY;
49			}
50
51			//* This section identical for most OBMol conversions *
52			////////////////////////////////////////////////////
53			/// The "API" interface functions
54			virtual bool WriteMolecule(OBBase* pOb, OBConversion* pConv);
55
56			private:
57			bool AreSameFragments(OBMol& mol, vector<int>& frag1, vector<int>& frag2);
58			void findAngles(OBMol& mol);
59			OBMol* createMolFromFragment(OBMol& mol, vector<int>& fragment);
60			void WriteMDFile(vector<OBMol*> mols, vector<int> numMols, ostream& os, OBMol& mol, vector<int>& indices);
61			};
62
63			//Make an instance of the format class
64			OOPSEFormat theOOPSEFormat;
65
66			/////////////////////////////////////////////////////////////////
67
68
69			bool OOPSEFormat::WriteMolecule(OBBase* pOb, OBConversion* pConv)
70			{
71			OBMol* pmol = dynamic_cast<OBMol*>(pOb);
72			if(pmol==NULL)
73			return false;
74
75			vector<vector<int> > fragmentLists;
76			pmol->ContigFragList(fragmentLists);
77			OBBitVec unused;
78			vector<bool> used(fragmentLists.size(), 0);
79			vector<vector<int> > molecules;
80			vector<int> indices;
81			for(int i =0; i < used.size(); ++i) {
82			if (used[i])
83			{
84			continue;
85			}
86			used[i] = true;
87			vector<int> sameMolTypes;
88			sameMolTypes.push_back(i);
89			indices.insert(indices.end(), fragmentLists[i].begin(),
90			fragmentLists[i].end());
91			for (int j = i + 1;j < used.size(); ++j)
92			{
93			if (used[j])
94			{
95			continue;
96			}
97
98			if (AreSameFragments(*pmol, fragmentLists[i], fragmentLists[j]))
99			{
100			sameMolTypes.push_back(j);
101			indices.insert(indices.end(), fragmentLists[j].begin(),
102			fragmentLists[j].end());
103			used[j]=true;
104			}
105			}
106			molecules.push_back(sameMolTypes);
107
108			}
109
110			//
111			vector<OBMol*> mdMols;
112			vector<int> numMols;
113			for(vector<vector<int> >::iterator i = molecules.begin();
114			i != molecules.end(); ++i)
115			{
116			mdMols.push_back(createMolFromFragment(*pmol,
117			fragmentLists[i->front()]));
118			numMols.push_back((*i).size());
119			}
120
121			string OutputFileName = pConv->GetInFilename();
122			unsigned int pos = OutputFileName.rfind(".");
123			if(pos==string::npos)
124			OutputFileName += ".md";
125			else
126			OutputFileName = OutputFileName.substr(0, pos) + ".md";
127			ofstream ofs(OutputFileName.c_str());
128			if(!ofs)
129			{
130			cerr << "Cannot write to " << OutputFileName <<endl;
131			return false;
132			}
133
134			WriteMDFile(mdMols, numMols, ofs, *pmol, indices);
135
136			for(vector<OBMol*>::iterator i = mdMols.begin(); i != mdMols.end(); ++i)
137			{
138			delete *i;
139			}
140
141			//
142
143			return(true);
144			}
145
146			bool OOPSEFormat::AreSameFragments(OBMol& mol, vector<int>& frag1,
147			vector<int>& frag2)
148			{
149			if (frag1.size() != frag2.size())
150			{
151			return false;
152			}
153
154			//exact graph matching is a NP complete problem
155			/** @todo using sparse matrix to store the connectivities*/
156			for (unsigned int i =0 ; i < frag1.size(); ++i)
157			{
158			OBAtom* atom1 = mol.GetAtom(frag1[i]);
159			OBAtom* atom2 = mol.GetAtom(frag2[i]);
160
161			if (atom1->GetAtomicNum() != atom2->GetAtomicNum())
162			{
163			return false;
164			}
165			}
166			return true;
167			}
168
169			struct SameAngle
170			{
171			bool operator()(const triple<OBAtom,OBAtom,OBAtom*> t1,
172			const triple<OBAtom,OBAtom,OBAtom*> t2) const
173			{
174			return (t1.second == t2.second) && ( (t1.first == t2.first && t1.third == t2.third) \|\| (t1.first == t2.third && t1.third == t2.first));
175			}
176			};
177
178			void OOPSEFormat::findAngles(OBMol& mol) {
179			/*
180			//if already has data return
181			if(mol.HasData(OBGenericDataType::AngleData))
182			return;
183
184			vector<OBEdgeBase*>::iterator bi1,bi2;
185			OBBond* bond;
186			OBAtom a,b,*c;
187
188			set<triple<OBAtom,OBAtom,OBAtom*>, SameAngle> uniqueAngles;
189			//loop through all bonds generating torsions
190			for(bond = mol.BeginBond(bi1);bond;bond = mol.NextBond(bi1))
191			{
192			b = bond->GetBeginAtom();
193			c = bond->GetEndAtom();
194			if(b->IsHydrogen())
195			continue;
196
197			for(a = b->BeginNbrAtom(bi2);a;a = b->NextNbrAtom(bi2))
198			{
199			if(a == c)
200			continue;
201
202			uniqueAngles.insert(triple<OBAtom,OBAtom,OBAtom*>(a, b, c));
203			}
204			}
205
206			//get new data and attach it to molecule
207			OBAngleData *angles = new OBAngleData;
208			mol.SetData(angles);
209			set<triple<OBAtom,OBAtom,OBAtom*>, SameAngle>::iterator i;
210
211			for (i = uniqueAngles.begin(); i != uniqueAngles.end(); ++i) {
212			OBAngle angle;
213			angle.SetAtoms(i->first, i->second, i->second);
214			angles->SetData(angle);
215			}
216			*/
217			}
218
219			OBMol* OOPSEFormat::createMolFromFragment(OBMol& mol, vector<int>& fragment) {
220
221			OBMol* newMol = new OBMol();
222			newMol->ReserveAtoms(fragment.size());
223			newMol->BeginModify();
224			for(vector<int>::iterator i = fragment.begin(); i != fragment.end(); ++i) {
225			OBAtom* newAtom = newMol->NewAtom();
226			newAtom = mol.GetAtom(*i);
227			}
228			newMol->EndModify();
229			newMol->ConnectTheDots();
230			findAngles(*newMol);
231			newMol->FindTorsions();
232			return newMol;
233			}
234
235			void OOPSEFormat::WriteMDFile(vector<OBMol*> mols, vector<int> numMols, ostream& os, OBMol& mol, vector<int>& indices) {
236			std::string indentLevel1(" ");
237			std::string indentLevel2(" ");
238			std::string molPrefix("MolName");
239			unsigned int i;
240			const int BUFFLEN = 1024;
241			char buffer[BUFFLEN];
242			string str, str1;
243
244
245			os << "<OOPSE version=4>" << endl;
246			os << " <MetaData>" << endl << endl;
247
248			for(i = 0; i < mols.size(); ++i) {
249			OBMol* pmol = mols[i];
250
251			pmol->ConnectTheDots();
252			pmol->PerceiveBondOrders();
253			//pmol->FindSSSR();
254			//pmol->SetAromaticPerceived();
255			//pmol->Kekulize();
256
257			map<OBAtom*, int> atomMap;
258			os << "molecule {\n";
259			sprintf(buffer, "%d", i);
260			os << indentLevel1 << "name = " << "\"" << molPrefix << buffer << "\"" << ";\n";
261
262			//atom
263			int ai = 0;
264			FOR_ATOMS_OF_MOL(atom, *pmol ) {
265			str = atom->GetType();
266			ttab.SetFromType("INT");
267			ttab.SetToType("INT");
268			ttab.Translate(str1,str);
269			os << indentLevel1 << "atom[" << ai << "] {\n";
270			// os << indentLevel2 << "type = " << "\"" << etab.GetSymbol(atom->GetAtomicNum()) << "\"" << ";\n";
271			os << indentLevel2 << "type = " << "\"" << str1 << "\"" << ";\n";
272			os << indentLevel1 << "}\n";
273			atomMap[&(*atom)] = ai++;
274			}
275			os << "\n";
276
277			//bond
278			FOR_BONDS_OF_MOL(bond, *pmol ) {
279			os << indentLevel1 << "bond {\n";
280			os << indentLevel2 << "members(" << atomMap[bond->GetBeginAtom()] << ", " << atomMap[bond->GetEndAtom()] << ");\n";
281			os << indentLevel1 << "}\n";
282			}
283			/*
284			//bend
285			OBGenericData* pGenericData = pmol->GetData(OBGenericDataType::AngleData);
286			OBAngleData* pAngleData = dynamic_cast<OBAngleData*>(pGenericData);
287			vector<OBAngle> angles = pAngleData->GetData();
288
289			os << indentLevel1 << "nBends = " << angles.size() << ";\n";
290			int bendIndex = 0;
291			for (vector<OBAngle>::iterator ti = angles.begin(); ti != angles.end(); ++ti)
292			{
293			triple<OBAtom, OBAtom, OBAtom*> bendAtoms = ti->getAtoms();
294			os << indentLevel1 << "bend[" << bendIndex++ << "] {\n";
295			os << indentLevel2 << "member(" << atomMap[bendAtoms.first] << ", " << atomMap[bendAtoms.second] << atomMap[bendAtoms.third] <<");\n";
296			os << indentLevel1 << "}\n";
297			}
298
299			//torsion
300			pGenericData = pmol->GetData(OBGenericDataType::TorsionData);
301			OBTorsionData* pTorsionData = dynamic_cast<OBTorsionData*>(pGenericData);
302			vector<OBTorsion> torsions = pTorsionData->GetData();
303			vector<quad<OBAtom,OBAtom,OBAtom,OBAtom> > torsionArray;
304			for (vector<OBTorsion>::iterator ti = torsions.begin(); ti != torsions.end(); ++ti)
305			{
306			vector<quad<OBAtom,OBAtom,OBAtom,OBAtom> > tmpTorsions = ti->getTorsions();
307			torsionArray.insert(torsionArray.end(), tmpTorsions.begin(), tmpTorsions.end());
308			}
309
310			os << indentLevel1 << "nTorsions = " << torsionArray.size() << ";\n";
311			int torsionIndex = 0;
312			for (vector<quad<OBAtom,OBAtom,OBAtom,OBAtom> >::iterator ti = torsionArray.begin(); ti != torsionArray.end(); ++ti)
313			{
314			os << indentLevel1 << "torsion[" << torsionIndex++ << "] {\n";
315			os << indentLevel2 << "member(" << atomMap[ti->first] << ", " << atomMap[ti->second] <<", " << atomMap[ti->third] <<", " << atomMap[ti->forth] << ");\n";
316			os << indentLevel1 << "}\n";
317			}
318			*/
319			os << "}" << endl;
320			os << endl;
321
322			}
323
324			os << endl;
325
326
327			for(i=0; i < mols.size(); ++i) {
328			os << "component{" << endl;
329			sprintf(buffer, "%d", i);
330			os << indentLevel1 << "type = " << molPrefix << buffer << ";" << endl;
331			os << indentLevel1 << "nMol = " << numMols[i]<< ";" << endl;
332			os << "}" << endl;
333			}
334
335			os << " </MetaData>" << endl;
336			os << " <Snapshot>" << endl;
337			os << " <FrameData>" << endl;
338
339			sprintf(buffer, " Time: %.10g", 0.0);
340
341			os << buffer << endl;
342
343			sprintf(buffer, " Hmat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}", 100.0, 0.0, 0.0, 0.0, 100.0, 0.0, 0.0, 0.0, 100.0);
344
345			os << buffer << endl;
346			os << " </FrameData>" << endl;
347			os << " <StuntDoubles>" << endl;
348
349			OBAtom *atom;
350
351			for(vector<int>::iterator i = indices.begin();i != indices.end(); ++i) {
352			atom = mol.GetAtom(*i);
353			sprintf(buffer, "%10d %7s %18.10g %18.10g %18.10g %13e %13e %13e", *i - 1, "pv", atom->GetX(), atom->GetY(), atom->GetZ(), 0.0, 0.0, 0.0);
354			os << buffer << endl;
355			}
356			os << " </StuntDoubles>" << endl;
357			os << " </Snapshot>" << endl;
358			os << "</OOPSE>" << endl;
359			}
360
361			} //namespace OpenBabel
362