src/openbabel/xyzformat.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
 
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation version 2 of the License.
 
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
***********************************************************************/

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

#include <sstream>

using namespace std;
namespace OpenBabel
{

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

    virtual const char* Description() //required
    {
      return
        "XYZ cartesian coordinates format\n \
       Read Options e.g. -as\n\
        s  Output single bonds only\n\
        b  Disable bonding entirely\n\n";
    };

    virtual const char* SpecificationURL()
    {return "http://openbabel.sourceforge.net/formats/xyz.shtml";}; //optional

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

    //*** 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);
  };
  //***

  //Make an instance of the format class
  XYZFormat theXYZFormat;

  /////////////////////////////////////////////////////////////////
  bool XYZFormat::ReadMolecule(OBBase* pOb, OBConversion* pConv)
  {
    OBMol* pmol = dynamic_cast<OBMol*>(pOb);
    if(pmol==NULL)
      return false;

    //Define some references so we can use the old parameter names
    istream &ifs = *pConv->GetInStream();
    OBMol &mol = *pmol;
    const char* title = pConv->GetTitle();

    char buffer[BUFF_SIZE];

    stringstream errorMsg;

    unsigned int natoms;        // [ejk] assumed natoms could not be -ve

    if (!ifs.getline(buffer,BUFF_SIZE))
      {
        obErrorLog.ThrowError(__func__,
                              "Problems reading an XYZ file: Cannot read the first line.", obWarning);
        return(false);
      }

    if (sscanf(buffer, "%d", &natoms) == 0 || !natoms)
      {
        obErrorLog.ThrowError(__func__,
                              "Problems reading an XYZ file: The first line must contain the number of atoms.", obWarning);
        return(false);
      }

    mol.ReserveAtoms(natoms);

    // The next line contains a title string for the molecule. Use this
    // as the title for the molecule if the line is not
    // empty. Otherwise, use the title given by the calling function.
    if (!ifs.getline(buffer,BUFF_SIZE))
      {
        obErrorLog.ThrowError(__func__,
                              "Problems reading an XYZ file: Could not read the second line (title/comments).", obWarning);
        return(false);
      }
    string readTitle(buffer);
    string::size_type location = readTitle.find("Energy");
    if (location != string::npos)
      readTitle.erase(location);
    Trim(readTitle);

    location = readTitle.find_first_not_of(" \t\n\r");
    if (location != string::npos)
      mol.SetTitle(readTitle);
    else
      mol.SetTitle(title);

    mol.BeginModify();

    // The next lines contain four items each, separated by white
    // spaces: the atom type, and the coordinates of the atom
    vector<string> vs;
    for (unsigned int i = 1; i <= natoms; i ++)
      {
        if (!ifs.getline(buffer,BUFF_SIZE))
          {
            errorMsg << "Problems reading an XYZ file: "
                     << "Could not read line #" << i+2 << ", file error." << endl
                     << " According to line one, there should be " << natoms 
                     << " atoms, and therefore " << natoms+2 << " lines in the file.";

            obErrorLog.ThrowError(__func__, errorMsg.str() , obWarning);
            return(false);
          }
        tokenize(vs,buffer);
        if (vs.size() != 4)
          {
            errorMsg << "Problems reading an XYZ file: "
                     << "Could not read line #" << i+2 << "." << endl
                     << "OpenBabel found the line '" << buffer << "'" << endl
                     << "According to the specifications, this line should contain exactly 4 entries, separated by white space." << endl
                     << "However, OpenBabel found " << vs.size() << " items.";

            obErrorLog.ThrowError(__func__, errorMsg.str() , obWarning);
            return(false);
          }

        // Atom Type: get the atomic number from the element table, using
        // the first entry in the currently read line. If the entry makes
        // sense, set the atomic number and leave the atomic type open
        // (the type is then later faulted in when atom->GetType() is
        // called). If the entry does not make sense to use, set the atom
        // type manually, assuming that the author of the xyz-file had
        // something "special" in mind.
        OBAtom *atom  = mol.NewAtom();
        int atomicNum = etab.GetAtomicNum(vs[0].c_str());
        atom->SetAtomicNum(atomicNum); //set atomic number, or '0' if the atom type is not recognized
        if (atomicNum == 0)
          atom->SetType(vs[0]);

        // Read the atom coordinates
        char *endptr;
        double x = strtod((char*)vs[1].c_str(),&endptr);
        if (endptr == (char*)vs[1].c_str())
          {
            errorMsg << "Problems reading an XYZ file: "
                     << "Could not read line #" << i+2 << "." << endl
                     << "OpenBabel found the line '" << buffer << "'" << endl
                     << "According to the specifications, this line should contain exactly 4 entries, separated by white space." << endl
                     << "OpenBabel could not interpret item #1 as a number.";

            obErrorLog.ThrowError(__func__, errorMsg.str() , obWarning);
            return(false);
          }
        double y = strtod((char*)vs[2].c_str(),&endptr);
        if (endptr == (char*)vs[2].c_str())
          {
            errorMsg << "Problems reading an XYZ file: "
                     << "Could not read line #" << i+2 << "." << endl
                     << "OpenBabel found the line '" << buffer << "'" << endl
                     << "According to the specifications, this line should contain exactly 4 entries, separated by white space." << endl
                     << "OpenBabel could not interpret item #2 as a number.";

            obErrorLog.ThrowError(__func__, errorMsg.str() , obWarning);
            return(false);
          }
        double z = strtod((char*)vs[3].c_str(),&endptr);
        if (endptr == (char*)vs[3].c_str())
          {
            errorMsg << "Problems reading an XYZ file: "
                     << "Could not read line #" << i+2 << "." << endl
                     << "OpenBabel found the line '" << buffer << "'" << endl
                     << "According to the specifications, this line should contain exactly 4 entries, separated by white space." << endl
                     << "OpenBabel could not interpret item #3 as a number.";

            obErrorLog.ThrowError(__func__, errorMsg.str() , obWarning);
            return(false);
          }
        atom->SetVector(x,y,z); //set coordinates
      }

    // clean out any remaining blank lines
    while(ifs.peek() != EOF && ifs.good() && 
          (ifs.peek() == '\n' || ifs.peek() == '\r'))
      ifs.getline(buffer,BUFF_SIZE);

    if (!pConv->IsOption("b",OBConversion::INOPTIONS))
      mol.ConnectTheDots();
    if (!pConv->IsOption("s",OBConversion::INOPTIONS) && !pConv->IsOption("b",OBConversion::INOPTIONS))
      mol.PerceiveBondOrders();

    mol.EndModify();

    return(true);
  }

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

  bool XYZFormat::WriteMolecule(OBBase* pOb, OBConversion* pConv)
  {
    OBMol* pmol = dynamic_cast<OBMol*>(pOb);
    if(pmol==NULL)
      return false;

    //Define some references so we can use the old parameter names
    ostream &ofs = *pConv->GetOutStream();
    OBMol &mol = *pmol;

    unsigned int i;
    char buffer[BUFF_SIZE];

    snprintf(buffer, BUFF_SIZE, "%d", mol.NumAtoms());
    ofs << buffer << endl;
    if (fabs(mol.GetEnergy()) > 1.0e-3) // nonzero energy field
      snprintf(buffer, BUFF_SIZE, "%s\tEnergy: %15.7f",
               mol.GetTitle(), mol.GetEnergy());
    else
      snprintf(buffer, BUFF_SIZE, "%s", mol.GetTitle());
    ofs << buffer << endl;

    OBAtom *atom;
    string str,str1;
    for(i = 1;i <= mol.NumAtoms(); i++)
      {
        atom = mol.GetAtom(i);
        snprintf(buffer, BUFF_SIZE, "%3s%15.5f%15.5f%15.5f",
                 etab.GetSymbol(atom->GetAtomicNum()),
                 atom->GetX(),
                 atom->GetY(),
                 atom->GetZ());
        ofs << buffer << endl;
      }

    return(true);
  }

} //namespace OpenBabel
Revision:	1081
Committed:	Thu Oct 19 20:49:05 2006 UTC (19 years ago) by gezelter
File size:	9402 byte(s)
Log Message:	updated OpenBabel to version 2.0.2
#	User	Rev	Content
1	tim	741	/**********************************************************************
2			Copyright (C) 2000 by OpenEye Scientific Software, Inc.
3	gezelter	1081	Some portions Copyright (C) 2001-2006 by Geoffrey R. Hutchison
4	tim	741	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	gezelter	1081	#include "mol.hpp"
17			#include "obconversion.hpp"
18			#include "obmolecformat.hpp"
19	tim	741
20	gezelter	1081	#include <sstream>
21	tim	741
22			using namespace std;
23			namespace OpenBabel
24			{
25
26	gezelter	1081	class XYZFormat : public OBMoleculeFormat
27			{
28			public:
29			//Register this format type ID
30			XYZFormat()
31			{
32			OBConversion::RegisterFormat("xyz", this, "chemical/x-xyz");
33			}
34
35			virtual const char* Description() //required
36			{
37			return
38			"XYZ cartesian coordinates format\n \
39			Read Options e.g. -as\n\
40			s Output single bonds only\n\
41			b Disable bonding entirely\n\n";
42			};
43
44			virtual const char* SpecificationURL()
45			{return "http://openbabel.sourceforge.net/formats/xyz.shtml";}; //optional
46
47			virtual const char* GetMIMEType()
48			{ return "chemical/x-xyz"; };
49
50			//* This section identical for most OBMol conversions *
51			////////////////////////////////////////////////////
52			/// The "API" interface functions
53			virtual bool ReadMolecule(OBBase* pOb, OBConversion* pConv);
54			virtual bool WriteMolecule(OBBase* pOb, OBConversion* pConv);
55			};
56			//***
57
58			//Make an instance of the format class
59			XYZFormat theXYZFormat;
60
61			/////////////////////////////////////////////////////////////////
62			bool XYZFormat::ReadMolecule(OBBase* pOb, OBConversion* pConv)
63			{
64	tim	741	OBMol* pmol = dynamic_cast<OBMol*>(pOb);
65			if(pmol==NULL)
66	gezelter	1081	return false;
67	tim	741
68			//Define some references so we can use the old parameter names
69			istream &ifs = *pConv->GetInStream();
70			OBMol &mol = *pmol;
71			const char* title = pConv->GetTitle();
72	gezelter	1081
73	tim	741	char buffer[BUFF_SIZE];
74
75			stringstream errorMsg;
76
77			unsigned int natoms; // [ejk] assumed natoms could not be -ve
78
79			if (!ifs.getline(buffer,BUFF_SIZE))
80			{
81	gezelter	1081	obErrorLog.ThrowError(__func__,
82			"Problems reading an XYZ file: Cannot read the first line.", obWarning);
83			return(false);
84	tim	741	}
85
86			if (sscanf(buffer, "%d", &natoms) == 0 \|\| !natoms)
87	gezelter	1081	{
88			obErrorLog.ThrowError(__func__,
89			"Problems reading an XYZ file: The first line must contain the number of atoms.", obWarning);
90			return(false);
91			}
92	tim	741
93			mol.ReserveAtoms(natoms);
94
95			// The next line contains a title string for the molecule. Use this
96			// as the title for the molecule if the line is not
97			// empty. Otherwise, use the title given by the calling function.
98			if (!ifs.getline(buffer,BUFF_SIZE))
99	gezelter	1081	{
100			obErrorLog.ThrowError(__func__,
101			"Problems reading an XYZ file: Could not read the second line (title/comments).", obWarning);
102	tim	741	return(false);
103	gezelter	1081	}
104			string readTitle(buffer);
105			string::size_type location = readTitle.find("Energy");
106			if (location != string::npos)
107			readTitle.erase(location);
108			Trim(readTitle);
109
110			location = readTitle.find_first_not_of(" \t\n\r");
111			if (location != string::npos)
112			mol.SetTitle(readTitle);
113	tim	741	else
114	gezelter	1081	mol.SetTitle(title);
115	tim	741
116			mol.BeginModify();
117
118			// The next lines contain four items each, separated by white
119			// spaces: the atom type, and the coordinates of the atom
120			vector<string> vs;
121			for (unsigned int i = 1; i <= natoms; i ++)
122	gezelter	1081	{
123	tim	741	if (!ifs.getline(buffer,BUFF_SIZE))
124	gezelter	1081	{
125			errorMsg << "Problems reading an XYZ file: "
126			<< "Could not read line #" << i+2 << ", file error." << endl
127			<< " According to line one, there should be " << natoms
128			<< " atoms, and therefore " << natoms+2 << " lines in the file.";
129	tim	741
130	gezelter	1081	obErrorLog.ThrowError(__func__, errorMsg.str() , obWarning);
131			return(false);
132			}
133	tim	741	tokenize(vs,buffer);
134	gezelter	1081	if (vs.size() != 4)
135			{
136			errorMsg << "Problems reading an XYZ file: "
137			<< "Could not read line #" << i+2 << "." << endl
138			<< "OpenBabel found the line '" << buffer << "'" << endl
139			<< "According to the specifications, this line should contain exactly 4 entries, separated by white space." << endl
140			<< "However, OpenBabel found " << vs.size() << " items.";
141	tim	741
142	gezelter	1081	obErrorLog.ThrowError(__func__, errorMsg.str() , obWarning);
143			return(false);
144			}
145	tim	741
146			// Atom Type: get the atomic number from the element table, using
147			// the first entry in the currently read line. If the entry makes
148			// sense, set the atomic number and leave the atomic type open
149			// (the type is then later faulted in when atom->GetType() is
150			// called). If the entry does not make sense to use, set the atom
151			// type manually, assuming that the author of the xyz-file had
152			// something "special" in mind.
153			OBAtom *atom = mol.NewAtom();
154			int atomicNum = etab.GetAtomicNum(vs[0].c_str());
155			atom->SetAtomicNum(atomicNum); //set atomic number, or '0' if the atom type is not recognized
156			if (atomicNum == 0)
157	gezelter	1081	atom->SetType(vs[0]);
158	tim	741
159			// Read the atom coordinates
160			char *endptr;
161			double x = strtod((char*)vs[1].c_str(),&endptr);
162			if (endptr == (char*)vs[1].c_str())
163	gezelter	1081	{
164			errorMsg << "Problems reading an XYZ file: "
165			<< "Could not read line #" << i+2 << "." << endl
166			<< "OpenBabel found the line '" << buffer << "'" << endl
167			<< "According to the specifications, this line should contain exactly 4 entries, separated by white space." << endl
168			<< "OpenBabel could not interpret item #1 as a number.";
169	tim	741
170	gezelter	1081	obErrorLog.ThrowError(__func__, errorMsg.str() , obWarning);
171			return(false);
172			}
173	tim	741	double y = strtod((char*)vs[2].c_str(),&endptr);
174			if (endptr == (char*)vs[2].c_str())
175	gezelter	1081	{
176			errorMsg << "Problems reading an XYZ file: "
177			<< "Could not read line #" << i+2 << "." << endl
178			<< "OpenBabel found the line '" << buffer << "'" << endl
179			<< "According to the specifications, this line should contain exactly 4 entries, separated by white space." << endl
180			<< "OpenBabel could not interpret item #2 as a number.";
181	tim	741
182	gezelter	1081	obErrorLog.ThrowError(__func__, errorMsg.str() , obWarning);
183			return(false);
184			}
185	tim	741	double z = strtod((char*)vs[3].c_str(),&endptr);
186			if (endptr == (char*)vs[3].c_str())
187	gezelter	1081	{
188			errorMsg << "Problems reading an XYZ file: "
189			<< "Could not read line #" << i+2 << "." << endl
190			<< "OpenBabel found the line '" << buffer << "'" << endl
191			<< "According to the specifications, this line should contain exactly 4 entries, separated by white space." << endl
192			<< "OpenBabel could not interpret item #3 as a number.";
193	tim	741
194	gezelter	1081	obErrorLog.ThrowError(__func__, errorMsg.str() , obWarning);
195			return(false);
196			}
197	tim	741	atom->SetVector(x,y,z); //set coordinates
198	gezelter	1081	}
199	tim	741
200			// clean out any remaining blank lines
201			while(ifs.peek() != EOF && ifs.good() &&
202	gezelter	1081	(ifs.peek() == '\n' \|\| ifs.peek() == '\r'))
203	tim	741	ifs.getline(buffer,BUFF_SIZE);
204
205			if (!pConv->IsOption("b",OBConversion::INOPTIONS))
206			mol.ConnectTheDots();
207			if (!pConv->IsOption("s",OBConversion::INOPTIONS) && !pConv->IsOption("b",OBConversion::INOPTIONS))
208			mol.PerceiveBondOrders();
209
210			mol.EndModify();
211
212			return(true);
213	gezelter	1081	}
214	tim	741
215	gezelter	1081	////////////////////////////////////////////////////////////////
216	tim	741
217	gezelter	1081	bool XYZFormat::WriteMolecule(OBBase* pOb, OBConversion* pConv)
218			{
219	tim	741	OBMol* pmol = dynamic_cast<OBMol*>(pOb);
220			if(pmol==NULL)
221	gezelter	1081	return false;
222	tim	741
223			//Define some references so we can use the old parameter names
224			ostream &ofs = *pConv->GetOutStream();
225			OBMol &mol = *pmol;
226
227			unsigned int i;
228			char buffer[BUFF_SIZE];
229
230	gezelter	1081	snprintf(buffer, BUFF_SIZE, "%d", mol.NumAtoms());
231	tim	741	ofs << buffer << endl;
232	gezelter	1081	if (fabs(mol.GetEnergy()) > 1.0e-3) // nonzero energy field
233			snprintf(buffer, BUFF_SIZE, "%s\tEnergy: %15.7f",
234			mol.GetTitle(), mol.GetEnergy());
235			else
236			snprintf(buffer, BUFF_SIZE, "%s", mol.GetTitle());
237	tim	741	ofs << buffer << endl;
238
239			OBAtom *atom;
240			string str,str1;
241			for(i = 1;i <= mol.NumAtoms(); i++)
242	gezelter	1081	{
243	tim	741	atom = mol.GetAtom(i);
244	gezelter	1081	snprintf(buffer, BUFF_SIZE, "%3s%15.5f%15.5f%15.5f",
245			etab.GetSymbol(atom->GetAtomicNum()),
246			atom->GetX(),
247			atom->GetY(),
248			atom->GetZ());
249	tim	741	ofs << buffer << endl;
250	gezelter	1081	}
251	tim	741
252			return(true);
253	gezelter	1081	}
254	tim	741
255			} //namespace OpenBabel