1 /**********************************************************************
2 Copyright (C) 2011 by Geoffrey Hutchison
3 
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation version 2 of the License.
7 
8 This program is distributed in the hope that it will be useful,
9 but WITHOUT ANY WARRANTY; without even the implied warranty of
10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
11 GNU General Public License for more details.
12 ***********************************************************************/
13 
14 #include <openbabel/babelconfig.h>
15 
16 #include <openbabel/obmolecformat.h>
17 #include <openbabel/mol.h>
18 #include <openbabel/atom.h>
19 #include <openbabel/bond.h>
20 #include <openbabel/obiter.h>
21 #include <openbabel/elements.h>
22 #include <openbabel/generic.h>
23 #include <cstdlib>
24 
25 
26 using namespace std;
27 namespace OpenBabel
28 {
29 
30   class XSFFormat : public OBMoleculeFormat
31   {
32   public:
33     //Register this format type ID
XSFFormat()34     XSFFormat()
35     {
36       OBConversion::RegisterFormat("xsf",this);
37       // animation variant
38       OBConversion::RegisterFormat("axsf",this);
39     }
40 
Description()41     virtual const char* Description() //required
42     {
43       return
44         "XCrySDen Structure Format\n"
45         "Read Options e.g. -as\n"
46         "  s  Output single bonds only\n"
47         "  b  Disable bonding entirely\n\n";
48     };
49 
SpecificationURL()50     virtual const char* SpecificationURL()
51     { return "http://www.xcrysden.org/doc/XSF.html/" ;}; //optional
52 
53     //Flags() can return be any the following combined by | or be omitted if none apply
54     // NOTREADABLE  READONEONLY  NOTWRITABLE  WRITEONEONLY
Flags()55     virtual unsigned int Flags()
56     {
57       return READONEONLY | NOTWRITABLE;
58     };
59 
60     //*** This section identical for most OBMol conversions ***
61     ////////////////////////////////////////////////////
62     /// The "API" interface functions
63     virtual bool ReadMolecule(OBBase* pOb, OBConversion* pConv);
64   };
65   //***
66 
67   //Make an instance of the format class
68   XSFFormat theXSFFormat;
69 
70   /////////////////////////////////////////////////////////////////
ReadMolecule(OBBase * pOb,OBConversion * pConv)71   bool XSFFormat::ReadMolecule(OBBase* pOb, OBConversion* pConv)
72   {
73 
74     OBMol* pmol = pOb->CastAndClear<OBMol>();
75     if (pmol == nullptr)
76       return false;
77 
78     //Define some references so we can use the old parameter names
79     istream &ifs = *pConv->GetInStream();
80     OBMol &mol = *pmol;
81     const char* title = pConv->GetTitle();
82 
83     char buffer[BUFF_SIZE];
84     string str;
85     double x,y,z;
86     OBAtom *atom;
87     vector3 translationVectors[3];
88     int numTranslationVectors = 0;
89     vector<string> vs;
90     vector<vector3> atomPositions;
91     bool createdAtoms = false;
92     int atomicNum;
93 
94     mol.BeginModify();
95 
96     while (ifs.getline(buffer, BUFF_SIZE))
97       {
98         if (buffer[0] == '#')
99           continue; // comment
100         if (strstr(buffer, "ATOMS") != nullptr) {
101           // Minimum of 4 columns -- AtNum, x, y, z (forces)
102           // where AtNum stands for atomic number (or symbol), while X Y Z are
103           ifs.getline(buffer, BUFF_SIZE);
104           tokenize(vs, buffer);
105           while (vs.size() >= 4) {
106             if (!createdAtoms) {
107               atom = mol.NewAtom();
108               //set atomic number
109               atomicNum = OBElements::GetAtomicNum(vs[0].c_str());
110               if (atomicNum == 0) {
111                 atomicNum = atoi(vs[0].c_str());
112               }
113               atom->SetAtomicNum(atomicNum);
114             }
115             x = atof((char*)vs[1].c_str());
116             y = atof((char*)vs[2].c_str());
117             z = atof((char*)vs[3].c_str());
118             atomPositions.push_back(vector3(x, y, z)); // we may have a movie or animation
119 
120             ifs.getline(buffer, BUFF_SIZE);
121             tokenize(vs, buffer);
122           }
123           createdAtoms = true; // don't run NewAtom() anymore
124         }
125         else if ( strstr(buffer, "PRIMVEC")
126                  || strstr(buffer, "CONVVEC") ) {
127           // translation vectors
128           numTranslationVectors = 0; // if we have an animation
129           while (numTranslationVectors < 3 && ifs.getline(buffer,BUFF_SIZE)) {
130             tokenize(vs,buffer); // we really need to check that it's 3 entries only
131             if (vs.size() < 3) return false; // timvdm 18/06/2008
132             x = atof((char*)vs[0].c_str());
133             y = atof((char*)vs[1].c_str());
134             z = atof((char*)vs[2].c_str());
135             translationVectors[numTranslationVectors++].Set(x, y, z);
136           }
137         }
138         else if (strstr(buffer, "PRIMCOORD") != nullptr) {
139           // read the coordinates
140           ifs.getline(buffer, BUFF_SIZE);
141           tokenize(vs, buffer);
142           if (vs.size() < 2) return false;
143           int numAtoms = atoi(vs[0].c_str());
144           for (int a = 0; a < numAtoms; ++a) {
145             if (!ifs.getline(buffer,BUFF_SIZE))
146               break;
147             tokenize(vs,buffer);
148             if (vs.size() < 4)
149               break;
150 
151             if (!createdAtoms) {
152               atom = mol.NewAtom();
153               //set atomic number
154               atomicNum = OBElements::GetAtomicNum(vs[0].c_str());
155               if (atomicNum == 0) {
156                 atomicNum = atoi(vs[0].c_str());
157               }
158               atom->SetAtomicNum(atomicNum);
159             }
160             x = atof((char*)vs[1].c_str());
161             y = atof((char*)vs[2].c_str());
162             z = atof((char*)vs[3].c_str());
163             atomPositions.push_back(vector3(x, y, z));
164           }
165         }
166       }
167 
168     mol.EndModify();
169 
170     int natom = mol.NumAtoms();
171     int numConformers = atomPositions.size() / natom;
172     for (int i = 0; i < numConformers; ++i) {
173       double *coordinates = new double[natom * 3];
174       for (int j = 0; j < natom; ++j) {
175         vector3 currentPosition = atomPositions[i*natom + j];
176         coordinates[j*3] = currentPosition.x();
177         coordinates[j*3 + 1] = currentPosition.y();
178         coordinates[j*3 + 2] = currentPosition.z();
179       }
180       mol.AddConformer(coordinates);
181     }
182     // Delete first conformer, created by EndModify, bunch of 0s
183     mol.DeleteConformer(0);
184     // Set geometry to last one
185     mol.SetConformer(mol.NumConformers() - 1);
186 
187     if (!pConv->IsOption("b",OBConversion::INOPTIONS))
188       mol.ConnectTheDots();
189     if (!pConv->IsOption("s",OBConversion::INOPTIONS) && !pConv->IsOption("b",OBConversion::INOPTIONS))
190       mol.PerceiveBondOrders();
191 
192     // Add final properties
193     mol.SetTitle(title);
194     if (numTranslationVectors == 3) {
195       OBUnitCell *uc = new OBUnitCell;
196       uc->SetOrigin(fileformatInput);
197       uc->SetData(translationVectors[0],
198                   translationVectors[1],
199                   translationVectors[2]);
200       mol.SetData(uc);
201     }
202 
203     return(true);
204   }
205 
206 } //namespace OpenBabel
207