1 /**********************************************************************
2 Copyright (C) 2008 by Geoffrey R. Hutchison
3 Some portions Copyright (C) 2004-2008 by Chris Morley
4
5 This program is free software; you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation version 2 of the License.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
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/internalcoord.h>
23 #include <openbabel/generic.h>
24
25 #include <cstdlib>
26
27 //Possible replacement for strcasestr. See end of file
28 const char *_strcasestr(const char *s, const char *pattern);
29
30 #ifndef strcasestr
31 #define strcasestr _strcasestr
32 #endif
33
34 using namespace std;
35 namespace OpenBabel
36 {
37
38 class GaussianZMatrixInputFormat : public OBMoleculeFormat
39 {
40 public:
41 //Register this format type ID
GaussianZMatrixInputFormat()42 GaussianZMatrixInputFormat()
43 {
44 OBConversion::RegisterFormat("gzmat",this, "chemical/x-gaussian-input");
45 }
46
Description()47 virtual const char* Description() //required
48 {
49 return
50 "Gaussian Z-Matrix Input\n"
51 "Read Options e.g. -as\n"
52 " s Output single bonds only\n"
53 " b Disable bonding entirely\n\n"
54 "Write Options e.g. -xk\n"
55 " k \"keywords\" Use the specified keywords for input\n"
56 " f <file> Read the file specified for input keywords\n\n";
57 };
58
SpecificationURL()59 virtual const char* SpecificationURL()
60 { return "https://www.gaussian.com/zmat/"; };
61
GetMIMEType()62 virtual const char* GetMIMEType()
63 { return "chemical/x-gaussian-input"; };
64
65 /// The "API" interface functions
66 virtual bool ReadMolecule(OBBase* pOb, OBConversion* pConv);
67 virtual bool WriteMolecule(OBBase* pOb, OBConversion* pConv);
68 };
69
70 //Make an instance of the format class
71 GaussianZMatrixInputFormat theGaussianZMatrixInputFormat;
72
73 ////////////////////////////////////////////////////////////////
74
WriteMolecule(OBBase * pOb,OBConversion * pConv)75 bool GaussianZMatrixInputFormat::WriteMolecule(OBBase* pOb, OBConversion* pConv)
76 {
77 OBMol* pmol = dynamic_cast<OBMol*>(pOb);
78 if (pmol == nullptr)
79 return false;
80
81 //Define some references so we can use the old parameter names
82 ostream &ofs = *pConv->GetOutStream();
83 OBMol &mol = *pmol;
84
85 char buffer[BUFF_SIZE];
86 const char *keywords = pConv->IsOption("k",OBConversion::OUTOPTIONS);
87 const char *keywordsEnable = pConv->IsOption("k",OBConversion::GENOPTIONS);
88 const char *keywordFile = pConv->IsOption("f",OBConversion::OUTOPTIONS);
89 string defaultKeywords = "!Put Keywords Here, check Charge and Multiplicity.\n#";
90
91 if(keywords) {
92 defaultKeywords = keywords;
93 }
94
95 if (keywordsEnable) {
96 string model;
97 string basis;
98 string method;
99
100 OBPairData *pd = (OBPairData *) pmol->GetData("model");
101 if(pd)
102 model = pd->GetValue();
103
104 pd = (OBPairData *) pmol->GetData("basis");
105 if(pd)
106 basis = pd->GetValue();
107
108 pd = (OBPairData *) pmol->GetData("method");
109 if(pd)
110 method = pd->GetValue();
111
112 if(method == "optimize") {
113 method = "opt";
114 }
115
116 if(model != "" && basis != "" && method != "") {
117 ofs << model << "/" << basis << "," << method << "\n";
118 }
119 else {
120 ofs << "#Unable to translate keywords!\n";
121 ofs << defaultKeywords << "\n";
122 }
123 }
124 else if (keywordFile) {
125 ifstream kfstream(keywordFile);
126 string keyBuffer;
127 if (kfstream) {
128 while (getline(kfstream, keyBuffer))
129 ofs << keyBuffer << "\n";
130 }
131 }
132 else {
133 ofs << defaultKeywords << "\n";
134 }
135 ofs << "\n"; // blank line after keywords
136 ofs << " " << mol.GetTitle() << "\n\n";
137
138 snprintf(buffer, BUFF_SIZE, "%d %d",
139 mol.GetTotalCharge(),
140 mol.GetTotalSpinMultiplicity());
141 ofs << buffer << "\n";
142
143 // OK, now that we set the keywords, charge, etc. we generate the z-matrix
144
145 OBAtom *a,*b,*c;//, *atom;
146
147 vector<OBInternalCoord*> vic;
148 vic.push_back(nullptr);
149 FOR_ATOMS_OF_MOL(atom, mol)
150 vic.push_back(new OBInternalCoord);
151
152 CartesianToInternal(vic,mol);
153
154 double r,w,t;
155 string type;
156 FOR_ATOMS_OF_MOL(atom, mol)
157 {
158 a = vic[atom->GetIdx()]->_a;
159 b = vic[atom->GetIdx()]->_b;
160 c = vic[atom->GetIdx()]->_c;
161
162 type = OBElements::GetSymbol(atom->GetAtomicNum());
163 if (atom->GetIsotope() != 0) {
164 snprintf(buffer, BUFF_SIZE, "(Iso=%d)", atom->GetIsotope());
165 type += buffer;
166 }
167
168 switch(atom->GetIdx())
169 {
170 case 1:
171 snprintf(buffer, BUFF_SIZE, "%-s\n",type.c_str());
172 ofs << buffer;
173 continue;
174 break;
175
176 case 2:
177 snprintf(buffer, BUFF_SIZE, "%-s %d r%d\n",
178 type.c_str(), a->GetIdx(), atom->GetIdx());
179 ofs << buffer;
180 continue;
181 break;
182
183 case 3:
184 snprintf(buffer, BUFF_SIZE, "%-s %d r%d %d a%d\n",
185 type.c_str(), a->GetIdx(), atom->GetIdx(), b->GetIdx(), atom->GetIdx());
186 ofs << buffer;
187 continue;
188 break;
189
190 default:
191 snprintf(buffer, BUFF_SIZE, "%-s %d r%d %d a%d %d d%d\n",
192 type.c_str(), a->GetIdx(), atom->GetIdx(),
193 b->GetIdx(), atom->GetIdx(), c->GetIdx(), atom->GetIdx());
194 ofs << buffer;
195 }
196 }
197
198 ofs << "Variables:\n";
199
200 FOR_ATOMS_OF_MOL(atom, mol)
201 {
202 r = vic[atom->GetIdx()]->_dst;
203 w = vic[atom->GetIdx()]->_ang;
204 if (w < 0.0)
205 w += 360.0;
206 t = vic[atom->GetIdx()]->_tor;
207 if (t < 0.0)
208 t += 360.0;
209
210 switch(atom->GetIdx())
211 {
212 case 1:
213 continue;
214 break;
215
216 case 2:
217 snprintf(buffer, BUFF_SIZE, "r2= %6.4f\n", r);
218 ofs << buffer;
219 continue;
220 break;
221
222 case 3:
223 snprintf(buffer, BUFF_SIZE, "r3= %6.4f\na3= %6.2f\n", r, w);
224 ofs << buffer;
225 continue;
226 break;
227
228 default:
229 snprintf(buffer, BUFF_SIZE, "r%d= %6.4f\na%d= %6.2f\nd%d= %6.2f\n",
230 atom->GetIdx(), r, atom->GetIdx(), w, atom->GetIdx(), t);
231 ofs << buffer;
232 }
233 }
234
235 // file should end with a blank line
236 ofs << "\n";
237 return(true);
238 }
239
240 // This is based on the Babel 1.6 code. It may not be ideal.
241 // If you have problems (or examples of input files that do not work), please contact
242 // the openbabel-discuss@lists.sourceforge.net mailing list and/or post a bug
ReadMolecule(OBBase * pOb,OBConversion * pConv)243 bool GaussianZMatrixInputFormat::ReadMolecule(OBBase* pOb, OBConversion* pConv)
244 {
245 OBMol* pmol = pOb->CastAndClear<OBMol>();
246 if (pmol == nullptr)
247 return false;
248
249 //Define some references so we can use the old parameter names
250 istream &ifs = *pConv->GetInStream();
251 OBMol &mol = *pmol;
252 const char* title = pConv->GetTitle();
253
254 char buffer[BUFF_SIZE];
255
256 OBAtom *atom;
257 vector<OBInternalCoord*> vic;
258 vic.push_back(nullptr); // OBMol indexed from 1 -- potential atom index problem
259
260 vector<string> vs;
261 int charge = 0;
262 unsigned int spin = 1;
263 unsigned int blankLines = 0;
264 bool readVariables = false; // are we reading variables yet?
265 bool foundVariables = false; // have we seen actual text for a variable line?
266 map<string, double> variables; // map from variable name to value
267 vector<string> atomLines; // save atom lines to re-parse after reading variables
268
269 mol.BeginModify();
270
271 //@todo: Read keywords from this line
272 while (ifs.getline(buffer,BUFF_SIZE))
273 if (strncmp(buffer,"#", 1) == 0) // begins with '#'
274 break;
275
276 while (ifs.getline(buffer,BUFF_SIZE)) { // blank line, title, blank, then we begin
277 if (strlen(buffer) == 0)
278 blankLines++;
279 if (blankLines == 2)
280 break;
281 }
282
283 ifs.getline(buffer,BUFF_SIZE); // Charge = # Multiplicty = #
284 tokenize(vs, buffer, " \t\n");
285 if (vs.size() == 2)
286 {
287 charge = atoi(vs[0].c_str());
288 spin = atoi(vs[1].c_str());
289 }
290
291 // We read through atom lines and cache them (into atomLines)
292 while (ifs.getline(buffer,BUFF_SIZE)) {
293 if (strlen(buffer) == 0) {
294 if (foundVariables)
295 break; // blank line after variable section
296 else {
297 readVariables = true;
298 continue;
299 }
300 }
301
302 if (strcasestr(buffer, "VARIABLE") != nullptr) {
303 readVariables = true;
304 continue;
305 }
306
307 if (readVariables) {
308 tokenize(vs, buffer, "= \t\n");
309 if (vs.size() >= 2) {
310 variables[vs[0]] = atof(vs[1].c_str());
311 foundVariables = true;
312 // cerr << "var: " << vs[0] << " " << vs[1] << endl;
313 }
314 }
315 else {
316 atomLines.push_back(buffer);
317 vic.push_back(new OBInternalCoord);
318 }
319 } // end while
320
321 char *endptr;
322 double temp;
323 if (atomLines.size() > 0) {
324 unsigned int i, j;
325 for (i = 0; i < atomLines.size(); ++i) {
326 j = i+1;
327 tokenize(vs, atomLines[i]);
328 atom = mol.NewAtom();
329 atom->SetAtomicNum(OBElements::GetAtomicNum(vs[0].c_str()));
330
331 if (j == 1) {
332 continue; // first atom, just create it
333 }
334
335 if (j >= 2) {
336 if (vs.size() < 3) {return false;}
337 vic[j]->_a = mol.GetAtom(atoi(vs[1].c_str()));
338
339 temp = strtod((char*)vs[2].c_str(), &endptr);
340 if (endptr != (char*)vs[2].c_str())
341 vic[j]->_dst = temp;
342 else
343 vic[j]->_dst = variables[vs[2].c_str()];
344 }
345
346 if (j >= 3) {
347 if (vs.size() < 5) {return false;}
348 vic[j]->_b = mol.GetAtom(atoi(vs[3].c_str()));
349
350 temp = strtod((char*)vs[4].c_str(), &endptr);
351 if (endptr != (char*)vs[4].c_str())
352 vic[j]->_ang = temp;
353 else
354 vic[j]->_ang = variables[vs[4].c_str()];
355 }
356
357 if (j >= 4) {
358 if (vs.size() < 7) {return false;}
359 vic[j]->_c = mol.GetAtom(atoi(vs[5].c_str()));
360
361 temp = strtod((char*)vs[6].c_str(), &endptr);
362 if (endptr != (char*)vs[6].c_str()) {
363 vic[j]->_tor = temp;
364 } else {
365 const char* tor_str = vs[6].c_str();
366 if (tor_str[0] == '-')
367 vic[j]->_tor = -1 * variables[tor_str+1];
368 else
369 vic[j]->_tor = variables[tor_str];
370 }
371 }
372 }
373 }
374
375 if (mol.NumAtoms() == 0) { // e.g., if we're at the end of a file PR#1737209
376 mol.EndModify();
377 return false;
378 }
379
380 InternalToCartesian(vic,mol);
381
382 if (!pConv->IsOption("b",OBConversion::INOPTIONS))
383 mol.ConnectTheDots();
384 if (!pConv->IsOption("s",OBConversion::INOPTIONS) && !pConv->IsOption("b",OBConversion::INOPTIONS))
385 mol.PerceiveBondOrders();
386
387 mol.EndModify();
388
389 mol.SetTotalCharge(charge);
390 mol.SetTotalSpinMultiplicity(spin);
391
392 mol.SetTitle(title);
393 return(true);
394 }
395
396 } //namespace OpenBabel
397
398 //*****************************************************************************
399 // Adapted from the ht://Dig package <http://www.htdig.org/>
400 // Copyright (c) 1999, 2000, 2001 The ht://Dig Group
401 // For copyright details, see the file COPYING in your distribution
402 // or the GNU General Public License version 2 or later
403 // <http://www.gnu.org/copyleft/gpl.html>
404
405 const char *
_strcasestr(const char * s,const char * pattern)406 _strcasestr(const char *s, const char *pattern)
407 {
408 int length = strlen(pattern);
409
410 while (*s)
411 {
412 if (strncasecmp(s, pattern, length) == 0)
413 return s;
414 s++;
415 }
416 return nullptr;
417 }
418
419