applications/atom2md/openmdformat.cpp

/**********************************************************************
Copyright (C) 2000 by OpenEye Scientific Software, Inc.
Some portions Copyright (C) 2001-2006 by Geoffrey R. Hutchison
Some portions Copyright (C) 2004 by Chris Morley
Some portions Copyright (C) 2008-2009 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 <openbabel/obiter.h>
#include <openbabel/mol.h>
#include <openbabel/chains.h>
#include <openbabel/data.h>
#include <fstream>

#include "utils/StringUtils.hpp"

using namespace std;
namespace OpenBabel
{
  
  class OpenMDFormat : public OBMoleculeFormat
  {
  public:
    //Register this format type ID
    OpenMDFormat() 
    {      
      OBConversion::RegisterFormat("md",this);
    }
    
    virtual const char* Description() //required
    {
      return
        "OpenMD combined meta-data / cartesian coordinates format\n\
        No comments yet\n";
    };
    
    virtual const char* SpecificationURL()
    {return "http://www.openmd.net";}; //optional
    
    virtual const char* GetMIMEType() 
    {return "chemical/x-md"; };
    
    virtual unsigned int Flags() 
    { 
      return NOTREADABLE | WRITEONEONLY;
    }
    
    virtual bool WriteMolecule(OBBase* pOb, OBConversion* pConv);

  private:
    bool AreSameFragments(OBMol& mol, vector<int>& frag1, vector<int>& frag2);
    OBMol* createMolFromFragment(OBMol& mol, vector<int>& fragment);
    void WriteMDFile(vector<OBMol*> mols, vector<int> numMols, ostream& os, 
                     OBMol& mol, vector<int>& indices);
    void CalcBoundingBox(OBMol &mol,
                         double &min_x, double &max_x,
                         double &min_y, double &max_y,
                         double &min_z, double &max_z);
      
  };
  
  //Make an instance of the format class
  OpenMDFormat theOpenMDFormat;
 
  bool OpenMDFormat::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();
    size_t pos = OutputFileName.rfind(".");
    if(pos!=string::npos)
      OutputFileName = OutputFileName.substr(0, pos) + ".md";       
    else
      OutputFileName += ".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 OpenMDFormat::AreSameFragments(OBMol& mol, vector<int>& frag1, 
                                     vector<int>& frag2) {
    if (frag1.size() != frag2.size())
      return false;
    
    // Exact graph matching is an NP complete problem.
    // This just matches all of the atom atomic numbers and may falsely 
    // detect identical fragments which aren't really identical.
    // @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));
    }
  };


  OBMol* OpenMDFormat::createMolFromFragment(OBMol& mol, 
                                            vector<int>& fragment) {
    
    OBMol* newMol = new OBMol();
    OBChainsParser* chainParser = new OBChainsParser();   
    bool molIsWater = false;

    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();
    newMol->PerceiveBondOrders();

    return newMol;
  }
  
  void OpenMDFormat::WriteMDFile(vector<OBMol*> mols, vector<int> numMols, 
                                ostream& os, OBMol& mol, 
                                vector<int>& indices) {
    
    std::string molPrefix("MolName");
    std::string resName;
    unsigned int i;
    const int BUFFLEN = 1024;
    char buffer[BUFFLEN];
    string str, str1, str2, str3;
    OBAtom *a, *b, *c, *d;    
    bool molIsWater = false;
    OBResidue *r;
    int resKey, myserial;
    char type_name[10];
    char *element_name;
    int res_num;
    OBChainsParser* chainParser = new OBChainsParser();   
    double min_x, max_x, min_y, max_y, min_z, max_z; /* Edges of bounding box */
    
    os << "<OpenMD version=1>" << endl;
    os << "  <MetaData>" << endl << endl;
    
    for(i = 0; i < mols.size(); ++i) {
      OBMol* pmol = mols[i];
      map<OBAtom*, int> atomMap;


      //chainParser->PerceiveChains(*pmol, false);
      molIsWater = false;
      FOR_RESIDUES_OF_MOL(residue, *pmol) {
        if (residue->GetName().compare("HOH") == 0) {
          molIsWater = true;
        }
      }
      
      if (molIsWater) {
        // water include files define all of the known water types
        os << "#include \"water.md\";\n";
        pmol->SetTitle("HOH");
      } else {

        os << "molecule {\n";
        sprintf(buffer, "%d", i);
        os << "  name = \"" << molPrefix << buffer << "\";\n";
        
        int ai = 0;
        FOR_ATOMS_OF_MOL(atom, *pmol ) {
          str = atom->GetType();
          r = atom->GetResidue();
          
          if (r == NULL) 
            resName = "NULL";
          else 
            resName = r->GetName();
          
          if (resName.compare("NULL") ==0 || 
              resName.compare("LIG") == 0 || 
              resName.compare("UNK") == 0) {
            // Either couldn't find a residue at all or couldn't find a
            // reasonable residue name to use.  We'll punt and use
            // OpenBabel's internal atom typing:
            ttab.SetFromType("INT");
            ttab.SetToType("INT");
            ttab.Translate(str1, str);
          } else {                        
            
            // If we know what residue we've got, the specific atom name can
            // be used to help specify partial charges. 

            //resdat.SetResName(resName);
            
            // atom type from residue: 
            str = r->GetAtomID(&*atom);
           
            // arginine has separate indices for chemically-identical
            // nitrogen atoms:
            if (resName.compare("ARG") == 0) {
              if (str.compare("NH1") == 0 || str.compare("NH2") == 0) {
                str = "NH";
              }
            }
            if (resName.compare("VAL") == 0) {
              if (str.compare("CG1") == 0 || str.compare("CG2") == 0) {
                str = "CG";
              }
            }
            if (resName.compare("LEU") == 0) {
              if (str.compare("CD1") == 0 || str.compare("CD2") == 0) {
                str = "CD";
              }
            }
            if (resName.compare("ASP") == 0) {
              if (str.compare("OD1") == 0 || str.compare("OD2") == 0) {
                str = "OD";
              }
            }
            if (resName.compare("GLU") == 0) {
              if (str.compare("OE1") == 0 || str.compare("OE2") == 0) {
                str = "OE";
              }
            }
            if (resName.compare("TYR") == 0) {
              if (str.compare("CD1") == 0 || str.compare("CD2") == 0) {
                str = "CD";
              }
              if (str.compare("CE1") == 0 || str.compare("CE2") == 0) {
                str = "CE";
              }
            }
            

            if ((&*atom)->IsHydrogen()) {
               FOR_NBORS_OF_ATOM(nbr, *atom) {
                 str2 = r->GetAtomID(&*nbr);
                 size_t startpos = str2.find_first_not_of(" ");
                 size_t endpos = str2.find_last_not_of(" ");
                 if ((endpos - startpos) < 1) {
                   // if the bonded atom type has only one character (i.e. N)
                   // then the hydrogen will be labeled "HN" to show what
                   // kind of proton it is:
                   str3 = str2;
                 } else {
                   if (str2.compare("OH") == 0) {
                      str3 = "O";
                   } else {
                     // When the bonded atom type is more specific, we drop
                     // the first character:  i.e. H bonded to OG1 is HG1 type:
                     str3 = str2.substr(startpos+1, endpos-startpos);
                   }
                 }
                str = "H" + str3;
               }
               // same problem with arginine NH atoms, but now for connected hydrogens
               if (resName.compare("ARG") == 0) {
                 if (str.compare("HH1") == 0 || str.compare("HH2") == 0) {
                   str = "HH";
                 }
               }
               if (resName.compare("VAL") == 0) {
                 if (str.compare("HG1") == 0 || str.compare("HG2") == 0) {
                   str = "HG";
                 }
               }
               if (resName.compare("LEU") == 0) {
                 if (str.compare("HD1") == 0 || str.compare("HD2") == 0) {
                   str = "HD";
                 }
               }
               if (resName.compare("TYR") == 0) {
                 if (str.compare("HD1") == 0 || str.compare("HD2") == 0) {
                   str = "HD";
                 }
                 if (str.compare("HE1") == 0 || str.compare("HE2") == 0) {
                   str = "HE";
                 }
               }

            }

            // atom type from residue table:
            //resdat.LookupType(str, str2, hyb);
            size_t startpos = str.find_first_not_of(" ");
            size_t endpos = str.find_last_not_of(" ");
            str = str.substr( startpos, endpos-startpos+1 );
            str1 = resName + "-" + str;
          }      
          os << "  atom[" << ai << "] { ";
          os << "type = " << "\"" << str1 << "\"" << "; ";
          os << "position( " << (&*atom)->GetX() << ", " << (&*atom)->GetY() << ", " << (&*atom)->GetZ() << ");";
          os << "}\n";
          atomMap[&(*atom)] = ai++;
        }        
        os << "\n";
        
        //bond
        
        int b1, b2;
        FOR_BONDS_OF_MOL(bond, *pmol ) {
          b1 = atomMap[bond->GetBeginAtom()];
          b2 = atomMap[bond->GetEndAtom()];
          
          os << "  bond { ";
          
          if (b1 < b2) 
            os << "members(" << b1 <<  ", " << b2 << "); ";
          else
            os << "members(" << b2 <<  ", " << b1 << "); ";
          
          os << "}" << endl;
        } 
        
        os << endl;
       
        os << "}" << endl;
        os << endl;
      }
    }
    
    os << endl;
        
    for(i=0; i < mols.size(); ++i) {
      OBMol* pmol = mols[i];      
      os << "component{" << endl;
      if (std::string(pmol->GetTitle()).compare("HOH") == 0) {
        os << "  type = " << "\"HOH\"" << "; // change to appropriate water model" << endl;
      } else {
        sprintf(buffer, "%d", i);
        os << "  type = " << molPrefix << buffer << ";" << endl;
      }
      os << "  nMol = " << numMols[i]<< ";" << endl;
      os << "}" << endl;
    }
    
    os << "  </MetaData>" << endl;
    os << "  <Snapshot>" << endl;
    os << "    <FrameData>" << endl;
    
    sprintf(buffer, "        Time: %.10g", 0.0);
    
    os << buffer << endl;

    CalcBoundingBox(mol, min_x, max_x, min_y, max_y, min_z, max_z);

    // still to do: should compute a bounding box here
    sprintf(buffer, "        Hmat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}", 
            max_x - min_x, 0.0, 0.0, 0.0, max_y - min_y, 0.0, 0.0, 0.0, max_z - min_z);
    
    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 << "</OpenMD>" << endl;
  }

  void OpenMDFormat::CalcBoundingBox(OBMol &mol,
                                    double &min_x, double &max_x,
                                    double &min_y, double &max_y,
                                    double &min_z, double &max_z
                                    )
  {
    /* ---- Init bounding-box variables ---- */
    min_x = (double) 0.0;
    max_x = (double) 0.0;
    min_y = (double) 0.0;
    max_y = (double) 0.0;
    min_z = (double) 0.0;
    max_z = (double) 0.0;
    
    /* ---- Check all atoms ---- */
    for(unsigned int i = 1; i <= mol.NumAtoms(); ++i)
      {
        
        /* ---- Get a pointer to ith atom ---- */
        OBAtom *atom = mol.GetAtom(i);
        
        /* ---- Check for minimal/maximal x-position ---- */
        if (atom -> GetX() < min_x)
          min_x = atom -> GetX();
        if (atom -> GetX() > max_x)
          max_x = atom -> GetX();
        
        /* ---- Check for minimal/maximal y-position ---- */
        if (atom -> GetY() < min_y)
          min_y = atom -> GetY();
        if (atom -> GetY() > max_y)
          max_y = atom -> GetY();
        
        /* ---- Check for minimal/maximal z-position ---- */
        if (atom -> GetZ() < min_z)
          min_z = atom -> GetZ();
        if (atom -> GetZ() > max_z)
          max_z = atom -> GetZ();
        
      }    
  }
} //namespace OpenBabel

Revision:	1442
Committed:	Mon May 10 17:28:26 2010 UTC (14 years, 11 months ago) by gezelter
Original Path:	*trunk/src/applications/atom2md/openmdformat.cpp*
File size:	15317 byte(s)
Log Message:	Adding property set to svn entries
#	User	Rev	Content
1	gezelter	1390	/**********************************************************************
2			Copyright (C) 2000 by OpenEye Scientific Software, Inc.
3	gezelter	1399	Some portions Copyright (C) 2001-2006 by Geoffrey R. Hutchison
4	gezelter	1390	Some portions Copyright (C) 2004 by Chris Morley
5	gezelter	1399	Some portions Copyright (C) 2008-2009 by J. Daniel Gezelter
6	gezelter	1390
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 <openbabel/obiter.h>
20			#include <openbabel/mol.h>
21			#include <openbabel/chains.h>
22			#include <openbabel/data.h>
23			#include <fstream>
24
25			#include "utils/StringUtils.hpp"
26
27			using namespace std;
28			namespace OpenBabel
29			{
30
31			class OpenMDFormat : public OBMoleculeFormat
32			{
33			public:
34			//Register this format type ID
35			OpenMDFormat()
36			{
37			OBConversion::RegisterFormat("md",this);
38			}
39
40			virtual const char* Description() //required
41			{
42			return
43	gezelter	1399	"OpenMD combined meta-data / cartesian coordinates format\n\
44			No comments yet\n";
45	gezelter	1390	};
46
47			virtual const char* SpecificationURL()
48			{return "http://www.openmd.net";}; //optional
49
50			virtual const char* GetMIMEType()
51			{return "chemical/x-md"; };
52
53			virtual unsigned int Flags()
54			{
55			return NOTREADABLE \| WRITEONEONLY;
56			}
57
58			virtual bool WriteMolecule(OBBase* pOb, OBConversion* pConv);
59
60			private:
61			bool AreSameFragments(OBMol& mol, vector<int>& frag1, vector<int>& frag2);
62			OBMol* createMolFromFragment(OBMol& mol, vector<int>& fragment);
63			void WriteMDFile(vector<OBMol*> mols, vector<int> numMols, ostream& os,
64			OBMol& mol, vector<int>& indices);
65			void CalcBoundingBox(OBMol &mol,
66			double &min_x, double &max_x,
67			double &min_y, double &max_y,
68			double &min_z, double &max_z);
69
70			};
71
72			//Make an instance of the format class
73			OpenMDFormat theOpenMDFormat;
74
75			bool OpenMDFormat::WriteMolecule(OBBase* pOb, OBConversion* pConv) {
76			OBMol* pmol = dynamic_cast<OBMol*>(pOb);
77			if(pmol==NULL)
78			return false;
79
80			vector<vector<int> > fragmentLists;
81			pmol->ContigFragList(fragmentLists);
82			OBBitVec unused;
83			vector<bool> used(fragmentLists.size(), 0);
84			vector<vector<int> > molecules;
85			vector<int> indices;
86			for(int i =0; i < used.size(); ++i) {
87			if (used[i])
88			continue;
89
90			used[i] = true;
91			vector<int> sameMolTypes;
92			sameMolTypes.push_back(i);
93			indices.insert(indices.end(), fragmentLists[i].begin(),
94			fragmentLists[i].end());
95			for (int j = i + 1;j < used.size(); ++j) {
96			if (used[j])
97			continue;
98
99			if (AreSameFragments(*pmol, fragmentLists[i], fragmentLists[j])) {
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			vector<OBMol*> mdMols;
110			vector<int> numMols;
111			for(vector<vector<int> >::iterator i = molecules.begin();
112			i != molecules.end(); ++i) {
113
114			mdMols.push_back(createMolFromFragment(*pmol,
115			fragmentLists[i->front()]));
116			numMols.push_back((*i).size());
117			}
118
119			string OutputFileName = pConv->GetInFilename();
120			size_t pos = OutputFileName.rfind(".");
121			if(pos!=string::npos)
122			OutputFileName = OutputFileName.substr(0, pos) + ".md";
123			else
124			OutputFileName += ".md";
125
126			ofstream ofs(OutputFileName.c_str());
127			if(!ofs) {
128			cerr << "Cannot write to " << OutputFileName <<endl;
129			return false;
130			}
131
132	gezelter	1442
133
134	gezelter	1390	WriteMDFile(mdMols, numMols, ofs, *pmol, indices);
135
136			for(vector<OBMol*>::iterator i = mdMols.begin(); i != mdMols.end(); ++i) {
137			delete *i;
138			}
139
140			return(true);
141			}
142
143			bool OpenMDFormat::AreSameFragments(OBMol& mol, vector<int>& frag1,
144			vector<int>& frag2) {
145			if (frag1.size() != frag2.size())
146			return false;
147
148			// Exact graph matching is an NP complete problem.
149			// This just matches all of the atom atomic numbers and may falsely
150			// detect identical fragments which aren't really identical.
151			// @todo using sparse matrix to store the connectivities
152
153			for (unsigned int i =0 ; i < frag1.size(); ++i) {
154			OBAtom* atom1 = mol.GetAtom(frag1[i]);
155			OBAtom* atom2 = mol.GetAtom(frag2[i]);
156
157			if (atom1->GetAtomicNum() != atom2->GetAtomicNum())
158			return false;
159
160			}
161			return true;
162			}
163
164			struct SameAngle {
165			bool operator()(const triple<OBAtom,OBAtom,OBAtom*> t1,
166			const triple<OBAtom,OBAtom,OBAtom*> t2) const {
167			return (t1.second == t2.second) && ( (t1.first == t2.first && t1.third == t2.third) \|\| (t1.first == t2.third && t1.third == t2.first));
168			}
169			};
170
171
172			OBMol* OpenMDFormat::createMolFromFragment(OBMol& mol,
173			vector<int>& fragment) {
174
175			OBMol* newMol = new OBMol();
176	gezelter	1442	OBChainsParser* chainParser = new OBChainsParser();
177			bool molIsWater = false;
178
179	gezelter	1390	newMol->ReserveAtoms(fragment.size());
180			newMol->BeginModify();
181			for(vector<int>::iterator i = fragment.begin(); i != fragment.end(); ++i) {
182			OBAtom* newAtom = newMol->NewAtom();
183			newAtom = mol.GetAtom(*i);
184			}
185
186			newMol->EndModify();
187			newMol->ConnectTheDots();
188			newMol->PerceiveBondOrders();
189
190			return newMol;
191			}
192
193			void OpenMDFormat::WriteMDFile(vector<OBMol*> mols, vector<int> numMols,
194			ostream& os, OBMol& mol,
195			vector<int>& indices) {
196
197			std::string molPrefix("MolName");
198			std::string resName;
199			unsigned int i;
200			const int BUFFLEN = 1024;
201			char buffer[BUFFLEN];
202			string str, str1, str2, str3;
203			OBAtom a, b, c, d;
204			bool molIsWater = false;
205			OBResidue *r;
206			int resKey, myserial;
207			char type_name[10];
208			char *element_name;
209			int res_num;
210			OBChainsParser* chainParser = new OBChainsParser();
211			double min_x, max_x, min_y, max_y, min_z, max_z; /* Edges of bounding box */
212
213			os << "<OpenMD version=1>" << endl;
214			os << " <MetaData>" << endl << endl;
215
216			for(i = 0; i < mols.size(); ++i) {
217			OBMol* pmol = mols[i];
218			map<OBAtom*, int> atomMap;
219
220	gezelter	1442
221	gezelter	1398	//chainParser->PerceiveChains(*pmol, false);
222	gezelter	1390	molIsWater = false;
223			FOR_RESIDUES_OF_MOL(residue, *pmol) {
224			if (residue->GetName().compare("HOH") == 0) {
225			molIsWater = true;
226			}
227			}
228
229			if (molIsWater) {
230			// water include files define all of the known water types
231			os << "#include \"water.md\";\n";
232			pmol->SetTitle("HOH");
233			} else {
234
235			os << "molecule {\n";
236			sprintf(buffer, "%d", i);
237			os << " name = \"" << molPrefix << buffer << "\";\n";
238
239			int ai = 0;
240			FOR_ATOMS_OF_MOL(atom, *pmol ) {
241			str = atom->GetType();
242			r = atom->GetResidue();
243
244			if (r == NULL)
245			resName = "NULL";
246			else
247			resName = r->GetName();
248
249			if (resName.compare("NULL") ==0 \|\|
250			resName.compare("LIG") == 0 \|\|
251			resName.compare("UNK") == 0) {
252			// Either couldn't find a residue at all or couldn't find a
253			// reasonable residue name to use. We'll punt and use
254			// OpenBabel's internal atom typing:
255			ttab.SetFromType("INT");
256			ttab.SetToType("INT");
257			ttab.Translate(str1, str);
258			} else {
259
260			// If we know what residue we've got, the specific atom name can
261			// be used to help specify partial charges.
262
263			//resdat.SetResName(resName);
264
265			// atom type from residue:
266			str = r->GetAtomID(&*atom);
267
268			// arginine has separate indices for chemically-identical
269			// nitrogen atoms:
270			if (resName.compare("ARG") == 0) {
271			if (str.compare("NH1") == 0 \|\| str.compare("NH2") == 0) {
272			str = "NH";
273			}
274			}
275			if (resName.compare("VAL") == 0) {
276			if (str.compare("CG1") == 0 \|\| str.compare("CG2") == 0) {
277			str = "CG";
278			}
279			}
280			if (resName.compare("LEU") == 0) {
281			if (str.compare("CD1") == 0 \|\| str.compare("CD2") == 0) {
282			str = "CD";
283			}
284			}
285			if (resName.compare("ASP") == 0) {
286			if (str.compare("OD1") == 0 \|\| str.compare("OD2") == 0) {
287			str = "OD";
288			}
289			}
290			if (resName.compare("GLU") == 0) {
291			if (str.compare("OE1") == 0 \|\| str.compare("OE2") == 0) {
292			str = "OE";
293			}
294			}
295			if (resName.compare("TYR") == 0) {
296			if (str.compare("CD1") == 0 \|\| str.compare("CD2") == 0) {
297			str = "CD";
298			}
299			if (str.compare("CE1") == 0 \|\| str.compare("CE2") == 0) {
300			str = "CE";
301			}
302			}
303
304
305			if ((&*atom)->IsHydrogen()) {
306			FOR_NBORS_OF_ATOM(nbr, *atom) {
307			str2 = r->GetAtomID(&*nbr);
308			size_t startpos = str2.find_first_not_of(" ");
309			size_t endpos = str2.find_last_not_of(" ");
310			if ((endpos - startpos) < 1) {
311			// if the bonded atom type has only one character (i.e. N)
312			// then the hydrogen will be labeled "HN" to show what
313			// kind of proton it is:
314			str3 = str2;
315			} else {
316			if (str2.compare("OH") == 0) {
317			str3 = "O";
318			} else {
319			// When the bonded atom type is more specific, we drop
320			// the first character: i.e. H bonded to OG1 is HG1 type:
321			str3 = str2.substr(startpos+1, endpos-startpos);
322			}
323			}
324			str = "H" + str3;
325			}
326			// same problem with arginine NH atoms, but now for connected hydrogens
327			if (resName.compare("ARG") == 0) {
328			if (str.compare("HH1") == 0 \|\| str.compare("HH2") == 0) {
329			str = "HH";
330			}
331			}
332			if (resName.compare("VAL") == 0) {
333			if (str.compare("HG1") == 0 \|\| str.compare("HG2") == 0) {
334			str = "HG";
335			}
336			}
337			if (resName.compare("LEU") == 0) {
338			if (str.compare("HD1") == 0 \|\| str.compare("HD2") == 0) {
339			str = "HD";
340			}
341			}
342			if (resName.compare("TYR") == 0) {
343			if (str.compare("HD1") == 0 \|\| str.compare("HD2") == 0) {
344			str = "HD";
345			}
346			if (str.compare("HE1") == 0 \|\| str.compare("HE2") == 0) {
347			str = "HE";
348			}
349			}
350
351			}
352
353			// atom type from residue table:
354			//resdat.LookupType(str, str2, hyb);
355			size_t startpos = str.find_first_not_of(" ");
356			size_t endpos = str.find_last_not_of(" ");
357			str = str.substr( startpos, endpos-startpos+1 );
358			str1 = resName + "-" + str;
359	gezelter	1442	}
360	gezelter	1390	os << " atom[" << ai << "] { ";
361			os << "type = " << "\"" << str1 << "\"" << "; ";
362	gezelter	1397	os << "position( " << (&atom)->GetX() << ", " << (&atom)->GetY() << ", " << (&*atom)->GetZ() << ");";
363	gezelter	1390	os << "}\n";
364			atomMap[&(*atom)] = ai++;
365			}
366			os << "\n";
367
368			//bond
369
370			int b1, b2;
371			FOR_BONDS_OF_MOL(bond, *pmol ) {
372			b1 = atomMap[bond->GetBeginAtom()];
373			b2 = atomMap[bond->GetEndAtom()];
374
375			os << " bond { ";
376
377			if (b1 < b2)
378			os << "members(" << b1 << ", " << b2 << "); ";
379			else
380			os << "members(" << b2 << ", " << b1 << "); ";
381
382			os << "}" << endl;
383			}
384
385			os << endl;
386
387			os << "}" << endl;
388			os << endl;
389			}
390			}
391
392			os << endl;
393
394			for(i=0; i < mols.size(); ++i) {
395			OBMol* pmol = mols[i];
396			os << "component{" << endl;
397			if (std::string(pmol->GetTitle()).compare("HOH") == 0) {
398	gezelter	1442	os << " type = " << "\"HOH\"" << "; // change to appropriate water model" << endl;
399	gezelter	1390	} else {
400			sprintf(buffer, "%d", i);
401			os << " type = " << molPrefix << buffer << ";" << endl;
402			}
403			os << " nMol = " << numMols[i]<< ";" << endl;
404			os << "}" << endl;
405			}
406
407			os << " </MetaData>" << endl;
408			os << " <Snapshot>" << endl;
409			os << " <FrameData>" << endl;
410
411			sprintf(buffer, " Time: %.10g", 0.0);
412
413			os << buffer << endl;
414
415			CalcBoundingBox(mol, min_x, max_x, min_y, max_y, min_z, max_z);
416
417			// still to do: should compute a bounding box here
418			sprintf(buffer, " Hmat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}",
419			max_x - min_x, 0.0, 0.0, 0.0, max_y - min_y, 0.0, 0.0, 0.0, max_z - min_z);
420
421			os << buffer << endl;
422			os << " </FrameData>" << endl;
423			os << " <StuntDoubles>" << endl;
424	gezelter	1442
425	gezelter	1390	OBAtom *atom;
426	gezelter	1442
427	gezelter	1390	for(vector<int>::iterator i = indices.begin();i != indices.end(); ++i) {
428	gezelter	1442
429	gezelter	1390	atom = mol.GetAtom(*i);
430			sprintf(buffer, "%10d %7s %18.10g %18.10g %18.10g %13e %13e %13e", *i - 1,
431			"pv", atom->GetX(), atom->GetY(), atom->GetZ(), 0.0, 0.0, 0.0);
432			os << buffer << endl;
433			}
434			os << " </StuntDoubles>" << endl;
435			os << " </Snapshot>" << endl;
436			os << "</OpenMD>" << endl;
437			}
438
439			void OpenMDFormat::CalcBoundingBox(OBMol &mol,
440			double &min_x, double &max_x,
441			double &min_y, double &max_y,
442			double &min_z, double &max_z
443			)
444			{
445			/* ---- Init bounding-box variables ---- */
446			min_x = (double) 0.0;
447			max_x = (double) 0.0;
448			min_y = (double) 0.0;
449			max_y = (double) 0.0;
450			min_z = (double) 0.0;
451			max_z = (double) 0.0;
452
453			/* ---- Check all atoms ---- */
454			for(unsigned int i = 1; i <= mol.NumAtoms(); ++i)
455			{
456
457			/* ---- Get a pointer to ith atom ---- */
458			OBAtom *atom = mol.GetAtom(i);
459
460			/* ---- Check for minimal/maximal x-position ---- */
461			if (atom -> GetX() < min_x)
462			min_x = atom -> GetX();
463			if (atom -> GetX() > max_x)
464			max_x = atom -> GetX();
465
466			/* ---- Check for minimal/maximal y-position ---- */
467			if (atom -> GetY() < min_y)
468			min_y = atom -> GetY();
469			if (atom -> GetY() > max_y)
470			max_y = atom -> GetY();
471
472			/* ---- Check for minimal/maximal z-position ---- */
473			if (atom -> GetZ() < min_z)
474			min_z = atom -> GetZ();
475			if (atom -> GetZ() > max_z)
476			max_z = atom -> GetZ();
477
478			}
479			}
480			} //namespace OpenBabel
481