1 /********************************************************************** 2 Copyright (C) 2000 by OpenEye Scientific Software, Inc. 3 Some portions Copyright (C) 2001-2010 by Geoffrey R. Hutchison 4 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 #include <openbabel/babelconfig.h> 16 17 #include <openbabel/data.h> 18 #include <openbabel/data_utilities.h> 19 #include <openbabel/obmolecformat.h> 20 #include <openbabel/mol.h> 21 #include <openbabel/atom.h> 22 #include <openbabel/bond.h> 23 #include <openbabel/obiter.h> 24 #include <openbabel/elements.h> 25 #include <openbabel/generic.h> 26 27 #include <openbabel/pointgroup.h> 28 #include <cstdlib> 29 30 using namespace std; 31 namespace OpenBabel 32 { 33 34 class GaussianOutputFormat : public OBMoleculeFormat 35 { 36 public: 37 //Register this format type ID GaussianOutputFormat()38 GaussianOutputFormat() 39 { 40 OBConversion::RegisterFormat("gal",this, "chemical/x-gaussian-log"); 41 OBConversion::RegisterFormat("g92",this); 42 OBConversion::RegisterFormat("g94",this); 43 OBConversion::RegisterFormat("g98",this); 44 OBConversion::RegisterFormat("g03",this); 45 OBConversion::RegisterFormat("g09",this); 46 OBConversion::RegisterFormat("g16",this); 47 } 48 Description()49 virtual const char* Description() //required 50 { 51 return 52 "Gaussian Output\n" 53 "Read Options e.g. -as\n" 54 " s Output single bonds only\n" 55 " b Disable bonding entirely\n\n"; 56 }; 57 SpecificationURL()58 virtual const char* SpecificationURL() 59 { return "https://www.gaussian.com/"; }; 60 GetMIMEType()61 virtual const char* GetMIMEType() 62 { return "chemical/x-gaussian-log"; }; 63 64 //Flags() can return be any the following combined by | or be omitted if none apply 65 // NOTREADABLE READONEONLY NOTWRITABLE WRITEONEONLY Flags()66 virtual unsigned int Flags() 67 { 68 return READONEONLY | NOTWRITABLE; 69 }; 70 71 /// The "API" interface functions 72 virtual bool ReadMolecule(OBBase* pOb, OBConversion* pConv); 73 }; 74 75 //Make an instance of the format class 76 GaussianOutputFormat theGaussianOutputFormat; 77 78 class GaussianInputFormat : public OBMoleculeFormat 79 { 80 public: 81 //Register this format type ID GaussianInputFormat()82 GaussianInputFormat() 83 { 84 OBConversion::RegisterFormat("com",this, "chemical/x-gaussian-input"); 85 OBConversion::RegisterFormat("gau",this); 86 OBConversion::RegisterFormat("gjc",this); 87 OBConversion::RegisterFormat("gjf",this); 88 OBConversion::RegisterOptionParam("b", nullptr, 0, OBConversion::OUTOPTIONS); 89 // Command-line keywords 90 OBConversion::RegisterOptionParam("k", nullptr, 1, OBConversion::OUTOPTIONS); 91 // Command-line keyword file 92 OBConversion::RegisterOptionParam("f", nullptr, 1, OBConversion::OUTOPTIONS); 93 } 94 Description()95 virtual const char* Description() //required 96 { 97 return 98 "Gaussian Input\n" 99 "Write Options e.g. -xk\n" 100 " b Output includes bonds\n" 101 " k \"keywords\" Use the specified keywords for input\n" 102 " f <file> Read the file specified for input keywords\n" 103 " u Write the crystallographic unit cell, if present.\n\n"; 104 }; 105 SpecificationURL()106 virtual const char* SpecificationURL() 107 { return "https://www.gaussian.com/input/"; }; 108 GetMIMEType()109 virtual const char* GetMIMEType() 110 { return "chemical/x-gaussian-input"; }; 111 112 //Flags() can return be any the following combined by | or be omitted if none apply 113 // NOTREADABLE READONEONLY NOTWRITABLE WRITEONEONLY Flags()114 virtual unsigned int Flags() 115 { 116 return NOTREADABLE | WRITEONEONLY; 117 }; 118 119 //////////////////////////////////////////////////// 120 /// The "API" interface functions 121 virtual bool WriteMolecule(OBBase* pOb, OBConversion* pConv); 122 123 }; 124 125 //Make an instance of the format class 126 GaussianInputFormat theGaussianInputFormat; 127 128 //////////////////////////////////////////////////////////////// 129 WriteMolecule(OBBase * pOb,OBConversion * pConv)130 bool GaussianInputFormat::WriteMolecule(OBBase* pOb, OBConversion* pConv) 131 { 132 OBMol* pmol = dynamic_cast<OBMol*>(pOb); 133 if (pmol == nullptr) 134 return false; 135 136 //Define some references so we can use the old parameter names 137 ostream &ofs = *pConv->GetOutStream(); 138 OBMol &mol = *pmol; 139 140 char buffer[BUFF_SIZE]; 141 const char *keywords = pConv->IsOption("k",OBConversion::OUTOPTIONS); 142 const char *keywordsEnable = pConv->IsOption("k",OBConversion::GENOPTIONS); 143 const char *keywordFile = pConv->IsOption("f",OBConversion::OUTOPTIONS); 144 bool writeUnitCell = (nullptr != pConv->IsOption("u", OBConversion::OUTOPTIONS)); 145 string defaultKeywords = "!Put Keywords Here, check Charge and Multiplicity.\n#"; 146 147 if(keywords) 148 { 149 defaultKeywords = keywords; 150 } 151 152 if (keywordsEnable) 153 { 154 string model; 155 string basis; 156 string method; 157 158 OBPairData *pd = (OBPairData *) pmol->GetData("model"); 159 if(pd) 160 model = pd->GetValue(); 161 162 pd = (OBPairData *) pmol->GetData("basis"); 163 if(pd) 164 basis = pd->GetValue(); 165 166 pd = (OBPairData *) pmol->GetData("method"); 167 if(pd) 168 method = pd->GetValue(); 169 170 if(method == "optimize") 171 { 172 method = "opt"; 173 } 174 175 if(model != "" && basis != "" && method != "") 176 { 177 ofs << model << "/" << basis << "," << method << endl; 178 } 179 else 180 { 181 ofs << "#Unable to translate keywords!" << endl; 182 ofs << defaultKeywords << endl; 183 } 184 } 185 else if (keywordFile) 186 { 187 ifstream kfstream(keywordFile); 188 string keyBuffer; 189 if (kfstream) 190 { 191 while (getline(kfstream, keyBuffer)) 192 ofs << keyBuffer << endl; 193 } 194 } 195 else 196 { 197 ofs << defaultKeywords << endl; 198 } 199 ofs << endl; // blank line after keywords 200 ofs << " " << mol.GetTitle() << endl << endl; 201 202 snprintf(buffer, BUFF_SIZE, "%d %d", 203 mol.GetTotalCharge(), 204 mol.GetTotalSpinMultiplicity()); 205 ofs << buffer << endl; 206 207 FOR_ATOMS_OF_MOL(atom, mol) 208 { 209 if (atom->GetIsotope() == 0) 210 snprintf(buffer, BUFF_SIZE, "%-3s %10.5f %10.5f %10.5f", 211 OBElements::GetSymbol(atom->GetAtomicNum()), 212 atom->GetX(), atom->GetY(), atom->GetZ()); 213 else 214 snprintf(buffer, BUFF_SIZE, "%-3s(Iso=%d) %10.5f %10.5f %10.5f", 215 OBElements::GetSymbol(atom->GetAtomicNum()), 216 atom->GetIsotope(), 217 atom->GetX(), atom->GetY(), atom->GetZ()); 218 219 ofs << buffer << endl; 220 } 221 // Translation vectors 222 OBUnitCell *uc = (OBUnitCell*)mol.GetData(OBGenericDataType::UnitCell); 223 if (uc && writeUnitCell) { 224 uc->FillUnitCell(&mol); // complete the unit cell with symmetry-derived atoms 225 226 vector<vector3> cellVectors = uc->GetCellVectors(); 227 for (vector<vector3>::iterator i = cellVectors.begin(); i != cellVectors.end(); ++i) { 228 snprintf(buffer, BUFF_SIZE, "TV %10.5f %10.5f %10.5f", 229 i->x(), 230 i->y(), 231 i->z()); 232 ofs << buffer << '\n'; 233 } 234 } 235 236 // Bonds, contributed by Daniel Mansfield 237 if (pConv->IsOption("b",OBConversion::OUTOPTIONS)) 238 { 239 // first, make begin.GetIdx < end.GetIdx 240 OBBond* bond; 241 OBAtom *atom; 242 vector<OBBond*>::iterator j; 243 vector<OBNodeBase*>::iterator i; 244 OBAtom *bgn, *end; 245 for (bond = mol.BeginBond(j); bond; bond = mol.NextBond(j)) 246 { 247 if (bond->GetBeginAtomIdx() > bond->GetEndAtomIdx()) { 248 bgn = bond->GetBeginAtom(); 249 end = bond->GetEndAtom(); 250 bond->SetBegin(end); 251 bond->SetEnd(bgn); 252 } 253 } 254 255 // this seems inefficient -- perhaps using atom neighbor iterators? 256 // -GRH 257 for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i)) 258 { 259 ofs << endl << atom->GetIdx() << " "; 260 for (bond = mol.BeginBond(j); bond; bond = mol.NextBond(j)) 261 { 262 if (bond->GetBeginAtomIdx() == atom->GetIdx()) { 263 snprintf(buffer, BUFF_SIZE, "%d %1.1f ", bond->GetEndAtomIdx(), (float) bond->GetBondOrder()); 264 ofs << buffer; 265 } 266 } 267 } // iterate through atoms 268 } // end writing bonds 269 270 // file should end with a blank line 271 ofs << endl; 272 return(true); 273 } 274 add_unique_pairdata_to_mol(OpenBabel::OBMol * mol,string attribute,string buffer,int start)275 static void add_unique_pairdata_to_mol(OpenBabel::OBMol *mol, 276 string attribute, 277 string buffer,int start) 278 { 279 int i; 280 vector<string> vs; 281 OpenBabel::OBPairData *pd; 282 string method; 283 284 tokenize(vs,buffer); 285 if (vs.size() >= start) 286 { 287 method = vs[start]; 288 for(i=start+1; (i<vs.size()); i++) 289 { 290 method.append(" "); 291 method.append(vs[i]); 292 } 293 pd = (OpenBabel::OBPairData *) mol->GetData(attribute); 294 if (nullptr == pd) 295 { 296 pd = new OpenBabel::OBPairData(); 297 pd->SetAttribute(attribute); 298 pd->SetOrigin(fileformatInput); 299 pd->SetValue(method); 300 mol->SetData(pd); 301 } 302 else 303 { 304 pd->SetValue(method); 305 } 306 } 307 } 308 extract_thermo(OpenBabel::OBMol * mol,string method,double temperature,double ezpe,double Hcorr,double Gcorr,double E0,double CV,int RotSymNum,std::vector<double> Scomponents)309 static int extract_thermo(OpenBabel::OBMol *mol,string method,double temperature, 310 double ezpe,double Hcorr,double Gcorr,double E0,double CV, 311 int RotSymNum,std::vector<double> Scomponents) 312 { 313 // Initiate correction database 314 OpenBabel::OBAtomicHeatOfFormationTable *ahof = new OpenBabel::OBAtomicHeatOfFormationTable(); 315 OpenBabel::OBAtomIterator OBai; 316 OpenBabel::OBAtom *OBa; 317 char valbuf[128]; 318 int ii,atomid,atomicnumber,found,foundall; 319 double dhofM0, dhofMT, S0MT, DeltaSMT; 320 double eFactor = HARTEE_TO_KCALPERMOL; 321 322 // Now loop over atoms in order to correct the Delta H formation 323 OBai = mol->BeginAtoms(); 324 atomid = 0; 325 foundall = 0; 326 dhofM0 = E0*eFactor; 327 dhofMT = dhofM0+(Hcorr-ezpe)*eFactor; 328 S0MT = 0; 329 if (temperature > 0) 330 { 331 // Multiply by 1000 to make the unit cal/mol K 332 S0MT += 1000*eFactor*(Hcorr-Gcorr)/temperature; 333 } 334 335 // Check for symmetry 336 OBPointGroup obPG; 337 338 obPG.Setup(mol); 339 const char *pg = obPG.IdentifyPointGroup(); 340 341 double Rgas = 1.9872041; // cal/mol K http://en.wikipedia.org/wiki/Gas_constant 342 double Srot = -Rgas * log(double(RotSymNum)); 343 344 345 //printf("DHf(M,0) = %g, DHf(M,T) = %g, S0(M,T) = %g\nPoint group = %s RotSymNum = %d Srot = %g\n", 346 // dhofM0, dhofMT, S0MT, pg, RotSymNum, Srot); 347 if (RotSymNum > 1) 348 { 349 // We assume Gaussian has done this correctly! 350 Srot = 0; 351 } 352 S0MT += Srot; 353 DeltaSMT = S0MT; 354 355 for (OBa = mol->BeginAtom(OBai); nullptr != OBa; OBa = mol->NextAtom(OBai)) 356 { 357 double dhfx0, dhfxT, S0xT; 358 atomicnumber = OBa->GetAtomicNum(); 359 found = ahof->GetHeatOfFormation(OBElements::GetSymbol(atomicnumber), 360 0, 361 method, 362 temperature, 363 &dhfx0, &dhfxT, &S0xT); 364 if (1 == found) 365 { 366 dhofM0 += dhfx0; 367 dhofMT += dhfxT; 368 DeltaSMT += S0xT; 369 foundall ++; 370 } 371 atomid++; 372 } 373 if (foundall == atomid) 374 { 375 std::string attr[5]; 376 double result[5]; 377 char buf[32]; 378 379 attr[0].assign("DeltaHform(0K)"); 380 result[0] = dhofM0; 381 snprintf(buf, sizeof(buf), "DeltaHform(%gK)", temperature); 382 attr[1].assign(buf); 383 result[1] = dhofMT; 384 snprintf(buf, sizeof(buf), "DeltaSform(%gK)", temperature); 385 attr[2].assign(buf); 386 result[2] = DeltaSMT; 387 snprintf(buf, sizeof(buf), "DeltaGform(%gK)", temperature); 388 attr[3].assign(buf); 389 result[3] = dhofMT - temperature*result[2]/1000; 390 snprintf(buf, sizeof(buf), "S0(%gK)", temperature); 391 attr[4].assign(buf); 392 result[4] = S0MT; 393 394 add_unique_pairdata_to_mol(mol, "method", method, 0); 395 for(ii=0; (ii<5); ii++) 396 { 397 // Add to molecule properties 398 sprintf(valbuf,"%f", result[ii]); 399 add_unique_pairdata_to_mol(mol, attr[ii], valbuf, 0); 400 } 401 sprintf(valbuf, "%f", ezpe*eFactor); 402 add_unique_pairdata_to_mol(mol, "zpe", valbuf, 0); 403 sprintf(valbuf, "%f", CV); 404 add_unique_pairdata_to_mol(mol, "cv", valbuf, 0); 405 sprintf(valbuf, "%f", CV+Rgas); 406 add_unique_pairdata_to_mol(mol, "cp", valbuf, 0); 407 // Entropy components 408 if (Scomponents.size() == 3) 409 { 410 const char *comps[3] = { "Strans", "Srot", "Svib" }; 411 for(int i=0; (i<3); i++) 412 { 413 sprintf(valbuf, "%f", Scomponents[i]); 414 add_unique_pairdata_to_mol(mol, comps[i], valbuf, 0); 415 } 416 } 417 // Finally store the energy in internal data structures as well. 418 mol->SetEnergy(dhofMT); 419 } 420 else 421 { 422 // Debug message? 423 } 424 // Clean up 425 delete ahof; 426 427 if (foundall == atomid) 428 return 1; 429 else 430 return 0; 431 } 432 433 // Reading Gaussian output has been tested for G98 and G03 to some degree 434 // If you have problems (or examples of older output), please contact 435 // the openbabel-discuss@lists.sourceforge.net mailing list and/or post a bug ReadMolecule(OBBase * pOb,OBConversion * pConv)436 bool GaussianOutputFormat::ReadMolecule(OBBase* pOb, OBConversion* pConv) 437 { 438 OBMol* pmol = pOb->CastAndClear<OBMol>(); 439 if (pmol == nullptr) 440 return false; 441 442 //Define some references so we can use the old parameter names 443 istream &ifs = *pConv->GetInStream(); 444 OBMol &mol = *pmol; 445 const char* title = pConv->GetTitle(); 446 447 char buffer[BUFF_SIZE]; 448 string str,str1,str2,thermo_method; 449 double x,y,z; 450 OBAtom *atom; 451 vector<string> vs,vs2; 452 int total_charge = 0; 453 unsigned int spin_multiplicity = 1; 454 bool hasPartialCharges = false; 455 string chargeModel; // descriptor for charges (e.g. "Mulliken") 456 457 // Variable for G2/G3/G4 etc. calculations 458 double ezpe,Hcorr,Gcorr,E0,CV; 459 bool ezpe_set=false,Hcorr_set=false,Gcorr_set=false,E0_set=false,CV_set=false; 460 double temperature = 0; /* Kelvin */ 461 std::vector<double> Scomponents; 462 // Electrostatic potential. ESP is calculated 463 // once unless the Opt and Pop jobs are combined. 464 // In this case, ESP is calculated once before 465 // the geometry optmization and once after. If this 466 // happens, the second ESP must be added to OBMol. 467 OBFreeGrid *esp = nullptr; 468 int NumEsp = 1; 469 int NumEspCounter = 0; 470 bool ESPisAdded = false; 471 472 // coordinates of all steps 473 // Set conformers to all coordinates we adopted 474 std::vector<double*> vconf; // index of all frames/conformers 475 std::vector<double> coordinates; // coordinates in each frame 476 int natoms = 0; // number of atoms -- ensure we don't go to a new job with a different molecule 477 478 // OBConformerData stores information about multiple steps 479 // we can change attribute later if needed (e.g., IRC) 480 OBConformerData *confData = new OBConformerData(); 481 confData->SetOrigin(fileformatInput); 482 std::vector<unsigned short> confDimensions = confData->GetDimension(); // to be fair, set these all to 3D 483 std::vector<double> confEnergies = confData->GetEnergies(); 484 std::vector< std::vector< vector3 > > confForces = confData->GetForces(); 485 486 //Vibrational data 487 std::vector< std::vector< vector3 > > Lx; 488 std::vector<double> Frequencies, Intensities; 489 //Rotational data 490 std::vector<double> RotConsts(3); 491 int RotSymNum=1; 492 OBRotationData::RType RotorType = OBRotationData::UNKNOWN; 493 494 // Translation vectors (if present) 495 vector3 translationVectors[3]; 496 int numTranslationVectors = 0; 497 498 //Electronic Excitation data 499 std::vector<double> Forces, Wavelengths, EDipole, 500 RotatoryStrengthsVelocity, RotatoryStrengthsLength; 501 502 // Orbital data 503 std::vector<double> orbitals; 504 std::vector<std::string> symmetries; 505 int aHOMO, bHOMO, betaStart; 506 aHOMO = bHOMO = betaStart = -1; 507 508 int i=0; 509 bool no_symmetry=false; 510 char coords_type[25]; 511 512 //Prescan file to find second instance of "orientation:" 513 //This will be the kind of coords used in the chk/fchk file 514 //Unless the "nosym" keyword has been requested 515 while (ifs.getline(buffer,BUFF_SIZE)) 516 { 517 if (strstr(buffer, "Symmetry turned off by external request.") != nullptr) 518 { 519 // The "nosym" keyword has been requested 520 no_symmetry = true; 521 } 522 if (strstr(buffer, "orientation:") != nullptr) 523 { 524 i++; 525 tokenize (vs, buffer); 526 // gotta check what types of orientation are present 527 strncpy (coords_type, vs[0].c_str(), 24); 528 strcat (coords_type, " orientation:"); 529 } 530 if ((no_symmetry && i==1) || i==2) 531 break; 532 } 533 // Reset end-of-file pointers etc. 534 ifs.clear(); 535 ifs.seekg(0); //rewind 536 537 mol.BeginModify(); 538 while (ifs.getline(buffer,BUFF_SIZE)) 539 { 540 if(strstr(buffer, "Entering Gaussian") != nullptr) 541 { 542 //Put some metadata into OBCommentData 543 string comment("Gaussian "); 544 545 if (nullptr != strchr(buffer, '=')) 546 { 547 comment += strchr(buffer,'=')+2; 548 comment += ""; 549 for(unsigned i=0; i<115 && ifs; ++i) 550 { 551 ifs.getline(buffer,BUFF_SIZE); 552 if (strstr(buffer, "Revision") != nullptr) 553 { 554 if (buffer[strlen(buffer)-1] == ',') 555 { 556 buffer[strlen(buffer)-1] = '\0'; 557 } 558 add_unique_pairdata_to_mol(&mol,"program",buffer,0); 559 } 560 else if(buffer[1]=='#') 561 { 562 //the line describing the method 563 if (strstr(buffer, "Opt") != nullptr) 564 { 565 // It is expected to have two sets of ESP in 566 // the log file if Opt is combined with Pop. 567 NumEsp = 2; 568 } 569 comment += buffer; 570 OBCommentData *cd = new OBCommentData; 571 cd->SetData(comment); 572 cd->SetOrigin(fileformatInput); 573 mol.SetData(cd); 574 575 tokenize(vs,buffer); 576 if (vs.size() > 1) 577 { 578 char *str = strdup(vs[1].c_str()); 579 char *ptr = strchr(str,'/'); 580 581 if (nullptr != ptr) 582 { 583 *ptr = ' '; 584 add_unique_pairdata_to_mol(&mol,"basis",ptr,0); 585 *ptr = '\0'; 586 add_unique_pairdata_to_mol(&mol,"method",str,0); 587 } 588 } 589 590 break; 591 } 592 } 593 } 594 } 595 596 else if (strstr(buffer, "Multiplicity") != nullptr) 597 { 598 tokenize(vs, buffer, " \t\n"); 599 if (vs.size() == 6) 600 { 601 total_charge = atoi(vs[2].c_str()); 602 spin_multiplicity = atoi(vs[5].c_str()); 603 } 604 605 ifs.getline(buffer,BUFF_SIZE); 606 } 607 else if (strstr(buffer, coords_type) != nullptr) 608 { 609 numTranslationVectors = 0; // ignore old translationVectors 610 ifs.getline(buffer,BUFF_SIZE); // --------------- 611 ifs.getline(buffer,BUFF_SIZE); // column headings 612 ifs.getline(buffer,BUFF_SIZE); // column headings 613 ifs.getline(buffer,BUFF_SIZE); // --------------- 614 ifs.getline(buffer,BUFF_SIZE); 615 tokenize(vs,buffer); 616 while (vs.size()>4) 617 { 618 int corr = vs.size()==5 ? -1 : 0; //g94; later versions have an extra column 619 x = atof((char*)vs[3+corr].c_str()); 620 y = atof((char*)vs[4+corr].c_str()); 621 z = atof((char*)vs[5+corr].c_str()); 622 int atomicNum = atoi((char*)vs[1].c_str()); 623 624 if (atomicNum > 0) // translation vectors are "-2" 625 { 626 if (natoms == 0) { // first time reading the molecule, create each atom 627 atom = mol.NewAtom(); 628 atom->SetAtomicNum(atoi((char*)vs[1].c_str())); 629 } 630 coordinates.push_back(x); 631 coordinates.push_back(y); 632 coordinates.push_back(z); 633 } 634 else { 635 translationVectors[numTranslationVectors++].Set(x, y, z); 636 } 637 638 if (!ifs.getline(buffer,BUFF_SIZE)) { 639 break; 640 } 641 tokenize(vs,buffer); 642 } 643 // done with reading atoms 644 natoms = mol.NumAtoms(); 645 if(natoms==0) 646 return false; 647 // malloc / memcpy 648 double *tmpCoords = new double [(natoms)*3]; 649 memcpy(tmpCoords, &coordinates[0], sizeof(double)*natoms*3); 650 vconf.push_back(tmpCoords); 651 coordinates.clear(); 652 confDimensions.push_back(3); // always 3D -- OBConformerData allows mixing 2D and 3D structures 653 } 654 else if(strstr(buffer, "Dipole moment") != nullptr) 655 { 656 ifs.getline(buffer,BUFF_SIZE); // actual components X ### Y #### Z ### 657 tokenize(vs,buffer); 658 if (vs.size() >= 6) 659 { 660 OBVectorData *dipoleMoment = new OBVectorData; 661 dipoleMoment->SetAttribute("Dipole Moment"); 662 double x, y, z; 663 x = atof(vs[1].c_str()); 664 y = atof(vs[3].c_str()); 665 z = atof(vs[5].c_str()); 666 dipoleMoment->SetData(x, y, z); 667 dipoleMoment->SetOrigin(fileformatInput); 668 mol.SetData(dipoleMoment); 669 } 670 if (!ifs.getline(buffer,BUFF_SIZE)) break; 671 } 672 else if (strstr(buffer, "Traceless Quadrupole moment") != nullptr) 673 { 674 ifs.getline(buffer,BUFF_SIZE); // actual components XX ### YY #### ZZ ### 675 tokenize(vs,buffer); 676 ifs.getline(buffer,BUFF_SIZE); // actual components XY ### XZ #### YZ ### 677 tokenize(vs2,buffer); 678 if ((vs.size() >= 6) && (vs2.size() >= 6)) 679 { 680 double Q[3][3]; 681 OpenBabel::OBMatrixData *quadrupoleMoment = new OpenBabel::OBMatrixData; 682 683 Q[0][0] = atof(vs[1].c_str()); 684 Q[1][1] = atof(vs[3].c_str()); 685 Q[2][2] = atof(vs[5].c_str()); 686 Q[1][0] = Q[0][1] = atof(vs2[1].c_str()); 687 Q[2][0] = Q[0][2] = atof(vs2[3].c_str()); 688 Q[2][1] = Q[1][2] = atof(vs2[5].c_str()); 689 matrix3x3 quad(Q); 690 691 quadrupoleMoment->SetAttribute("Traceless Quadrupole Moment"); 692 quadrupoleMoment->SetData(quad); 693 quadrupoleMoment->SetOrigin(fileformatInput); 694 mol.SetData(quadrupoleMoment); 695 } 696 if (!ifs.getline(buffer,BUFF_SIZE)) break; 697 } 698 else if (strstr(buffer, "Exact polarizability") != nullptr) 699 { 700 // actual components XX, YX, YY, XZ, YZ, ZZ 701 double xx, xy, yy, xz, yz, zz; 702 const char *ptr = buffer+strlen("Exact polarizability: "); 703 if (ptr && 704 6 == sscanf(ptr, "%8lf%8lf%8lf%8lf%8lf%8lf", 705 &xx, &xy, &yy, &xz, &yz, &zz)) 706 { 707 double Q[3][3]; 708 OpenBabel::OBMatrixData *pol_tensor = new OpenBabel::OBMatrixData; 709 710 Q[0][0] = xx; 711 Q[1][1] = yy; 712 Q[2][2] = zz; 713 Q[1][0] = Q[0][1] = xy; 714 Q[2][0] = Q[0][2] = xz; 715 Q[2][1] = Q[1][2] = yz; 716 matrix3x3 pol(Q); 717 718 if (mol.HasData("Exact polarizability")) 719 { 720 mol.DeleteData("Exact polarizability"); // Delete the old one to add the new one 721 } 722 pol_tensor->SetAttribute("Exact polarizability"); 723 pol_tensor->SetData(pol); 724 pol_tensor->SetOrigin(fileformatInput); 725 mol.SetData(pol_tensor); 726 } 727 if (!ifs.getline(buffer,BUFF_SIZE)) break; 728 } 729 else if(strstr(buffer, "Total atomic charges") != nullptr || 730 strstr(buffer, "Mulliken atomic charges") != nullptr || 731 strstr(buffer, "Mulliken charges:") != nullptr) 732 { 733 hasPartialCharges = true; 734 chargeModel = "Mulliken"; 735 /* 736 Gaussian usually calculates the electronic 737 properties more than once, before and after 738 geometry optimization. The second one is what 739 we should be interested in. Thus, here, we 740 delete the previously added Data to store the 741 new one. 742 */ 743 if (mol.HasData("Mulliken charges")) 744 { 745 mol.DeleteData("Mulliken charges"); 746 } 747 OBPcharge *Mulliken = new OpenBabel::OBPcharge(); 748 std::vector<double> MPA_q; 749 750 ifs.getline(buffer,BUFF_SIZE); // column headings 751 ifs.getline(buffer,BUFF_SIZE); 752 tokenize(vs,buffer); 753 while (vs.size() >= 3 && 754 strstr(buffer, "Sum of ") == nullptr) 755 { 756 atom = mol.GetAtom(atoi(vs[0].c_str())); 757 if (!atom) 758 break; 759 atom->SetPartialCharge(atof(vs[2].c_str())); 760 MPA_q.push_back(atof(vs[2].c_str())); 761 if (!ifs.getline(buffer,BUFF_SIZE)) break; 762 tokenize(vs,buffer); 763 764 } 765 if (MPA_q.size() == mol.NumAtoms()) 766 { 767 Mulliken->AddPartialCharge(MPA_q); 768 Mulliken->SetAttribute("Mulliken charges"); 769 Mulliken->SetOrigin(fileformatInput); 770 mol.SetData(Mulliken); 771 } 772 else 773 { 774 cout << "Read " << MPA_q.size() << " Mulliken charges for " << mol.NumAtoms() << " atoms\n"; 775 } 776 } 777 else if (strstr(buffer, "Hirshfeld charges") != nullptr && 778 strstr(buffer, "CM5 charges") != nullptr) 779 { 780 /* 781 Hirshfeld and CM5 charges are printed in the 782 same block in the Gaussian log file. 783 */ 784 hasPartialCharges = true; 785 chargeModel = "Hirshfeld"; 786 if (mol.HasData("Hirshfeld charges")) 787 { 788 mol.DeleteData("Hirshfeld charges"); 789 } 790 if (mol.HasData("CM5 charges")) 791 { 792 mol.DeleteData("CM5 charges"); 793 } 794 OBPcharge *Hirshfeld = new OpenBabel::OBPcharge(); 795 OBPcharge *CM5 = new OpenBabel::OBPcharge(); 796 std::vector<double> HPA_q; 797 std::vector<double> CM5_q; 798 ifs.getline(buffer,BUFF_SIZE); // column headings 799 ifs.getline(buffer,BUFF_SIZE); 800 tokenize(vs,buffer); 801 while (vs.size() >= 8 && 802 strstr(buffer, "Tot ") == nullptr) 803 { 804 atom = mol.GetAtom(atoi(vs[0].c_str())); 805 if (!atom) 806 break; 807 atom->SetPartialCharge(atof(vs[2].c_str())); 808 HPA_q.push_back(atof(vs[2].c_str())); 809 CM5_q.push_back(atof(vs[7].c_str())); 810 if (!ifs.getline(buffer,BUFF_SIZE)) break; 811 tokenize(vs,buffer); 812 813 } 814 if (CM5_q.size() == mol.NumAtoms() and 815 HPA_q.size() == mol.NumAtoms()) 816 { 817 Hirshfeld->AddPartialCharge(HPA_q); 818 Hirshfeld->SetAttribute("Hirshfeld charges"); 819 Hirshfeld->SetOrigin(fileformatInput); 820 CM5->AddPartialCharge(CM5_q); 821 CM5->SetAttribute("CM5 charges"); 822 CM5->SetOrigin(fileformatInput); 823 mol.SetData(CM5); 824 mol.SetData(Hirshfeld); 825 } 826 else 827 { 828 cout << "Read " << HPA_q.size() << " Hirshfeld charges for " << mol.NumAtoms() << " atoms\n"; 829 } 830 } 831 else if (strstr(buffer, "Electrostatic Properties Using The SCF Density") != nullptr) 832 { 833 NumEspCounter++; 834 } 835 else if (strstr(buffer, "Atomic Center") != nullptr && !ESPisAdded) 836 { 837 // Data points for ESP calculation 838 tokenize(vs,buffer); 839 if (nullptr == esp) 840 esp = new OpenBabel::OBFreeGrid(); 841 if (vs.size() == 8) 842 { 843 esp->AddPoint(atof(vs[5].c_str()),atof(vs[6].c_str()), 844 atof(vs[7].c_str()),0); 845 } 846 else if (vs.size() > 5) 847 { 848 double x,y,z; 849 if (3 == sscanf(buffer+32,"%10lf%10lf%10lf",&x,&y,&z)) 850 { 851 esp->AddPoint(x,y,z,0); 852 } 853 } 854 } 855 else if (strstr(buffer, "ESP Fit Center") != nullptr && !ESPisAdded) 856 { 857 // Data points for ESP calculation 858 tokenize(vs,buffer); 859 if (nullptr == esp) 860 esp = new OpenBabel::OBFreeGrid(); 861 if (vs.size() == 9) 862 { 863 esp->AddPoint(atof(vs[6].c_str()),atof(vs[7].c_str()), 864 atof(vs[8].c_str()),0); 865 } 866 else if (vs.size() > 6) 867 { 868 double x,y,z; 869 if (3 == sscanf(buffer+32,"%10lf%10lf%10lf",&x,&y,&z)) 870 { 871 esp->AddPoint(x,y,z,0); 872 } 873 } 874 } 875 else if (strstr(buffer, "Electrostatic Properties (Atomic Units)") != nullptr && !ESPisAdded) 876 { 877 int i,np; 878 OpenBabel::OBFreeGridPoint *fgp; 879 OpenBabel::OBFreeGridPointIterator fgpi; 880 for(i=0; (i<5); i++) 881 { 882 ifs.getline(buffer,BUFF_SIZE); // skip line 883 } 884 // Assume file is correct and that potentials are present 885 // where they should. 886 np = esp->NumPoints(); 887 fgpi = esp->BeginPoints(); 888 i = 0; 889 for(fgp = esp->BeginPoint(fgpi); nullptr != fgp; fgp = esp->NextPoint(fgpi)) 890 { 891 ifs.getline(buffer,BUFF_SIZE); 892 tokenize(vs,buffer); 893 if (vs.size() >= 2) 894 { 895 fgp->SetV(atof(vs[2].c_str())); 896 i++; 897 } 898 } 899 if (NumEsp == NumEspCounter) 900 { 901 if (i == np) 902 { 903 esp->SetAttribute("Electrostatic Potential"); 904 esp->SetOrigin(fileformatInput); 905 mol.SetData(esp); 906 ESPisAdded = true; 907 } 908 else 909 { 910 cout << "Read " << esp->NumPoints() << " ESP points i = " << i << "\n"; 911 } 912 } 913 else if (!ESPisAdded) 914 { 915 esp->Clear(); 916 } 917 } 918 else if (strstr(buffer, "Charges from ESP fit") != nullptr) 919 { 920 hasPartialCharges = true; 921 chargeModel = "ESP"; 922 if (mol.HasData("ESP charges")) 923 { 924 mol.DeleteData("ESP charges"); 925 } 926 OBPcharge *ESP = new OpenBabel::OBPcharge(); 927 std::vector<double> ESP_q; 928 ifs.getline(buffer,BUFF_SIZE); // Charge / dipole line 929 ifs.getline(buffer,BUFF_SIZE); // column header 930 ifs.getline(buffer,BUFF_SIZE); // real charges 931 tokenize(vs,buffer); 932 while (vs.size() >= 3 && 933 strstr(buffer, "-----") == nullptr) 934 { 935 atom = mol.GetAtom(atoi(vs[0].c_str())); 936 if (!atom) 937 break; 938 atom->SetPartialCharge(atof(vs[2].c_str())); 939 ESP_q.push_back(atof(vs[2].c_str())); 940 if (!ifs.getline(buffer,BUFF_SIZE)) break; 941 tokenize(vs,buffer); 942 } 943 if (ESP_q.size() == mol.NumAtoms()) 944 { 945 ESP->AddPartialCharge(ESP_q); 946 ESP->SetAttribute("ESP charges"); 947 ESP->SetOrigin(fileformatInput); 948 mol.SetData(ESP); 949 } 950 else 951 { 952 cout << "Read " << ESP_q.size() << " ESP charges for " << mol.NumAtoms() << " atoms\n"; 953 } 954 } 955 else if (strstr(buffer, "Natural Population") != nullptr) 956 { 957 hasPartialCharges = true; 958 chargeModel = "NBO"; 959 ifs.getline(buffer,BUFF_SIZE); // column headings 960 ifs.getline(buffer,BUFF_SIZE); // again 961 ifs.getline(buffer,BUFF_SIZE); // again (-----) 962 ifs.getline(buffer,BUFF_SIZE); // real data 963 tokenize(vs,buffer); 964 while (vs.size() >= 3 && 965 strstr(buffer, "=====") == nullptr) 966 { 967 atom = mol.GetAtom(atoi(vs[1].c_str())); 968 if (!atom) 969 break; 970 atom->SetPartialCharge(atof(vs[2].c_str())); 971 972 if (!ifs.getline(buffer,BUFF_SIZE)) break; 973 tokenize(vs,buffer); 974 } 975 } 976 else if(strstr(buffer, " Frequencies -- ")) //vibrational frequencies 977 { 978 //The info should appear only once as several blocks starting with this line 979 tokenize(vs, buffer); 980 for(unsigned int i=2; i<vs.size(); ++i) 981 Frequencies.push_back(atof(vs[i].c_str())); 982 ifs.getline(buffer,BUFF_SIZE); //Red. masses 983 ifs.getline(buffer,BUFF_SIZE); //Frc consts 984 ifs.getline(buffer,BUFF_SIZE); //IR Inten 985 tokenize(vs, buffer); 986 for(unsigned int i=3; i<vs.size(); ++i) 987 Intensities.push_back(atof(vs[i].c_str())); 988 989 ifs.getline(buffer, BUFF_SIZE); // column labels or Raman intensity 990 if(strstr(buffer, "Raman Activ")) { 991 ifs.getline(buffer, BUFF_SIZE); // Depolar (P) 992 ifs.getline(buffer, BUFF_SIZE); // Depolar (U) 993 ifs.getline(buffer, BUFF_SIZE); // column labels 994 } 995 ifs.getline(buffer, BUFF_SIZE); // actual displacement data 996 tokenize(vs, buffer); 997 vector<vector3> vib1, vib2, vib3; 998 double x, y, z; 999 while(vs.size() >= 5) { 1000 for (unsigned int i = 2; i < vs.size()-2; i += 3) { 1001 x = atof(vs[i].c_str()); 1002 y = atof(vs[i+1].c_str()); 1003 z = atof(vs[i+2].c_str()); 1004 1005 if (i == 2) 1006 vib1.push_back(vector3(x, y, z)); 1007 else if (i == 5) 1008 vib2.push_back(vector3(x, y, z)); 1009 else if (i == 8) 1010 vib3.push_back(vector3(x, y, z)); 1011 } 1012 1013 if (!ifs.getline(buffer, BUFF_SIZE)) 1014 break; 1015 tokenize(vs,buffer); 1016 } 1017 Lx.push_back(vib1); 1018 if (vib2.size()) 1019 Lx.push_back(vib2); 1020 if (vib3.size()) 1021 Lx.push_back(vib3); 1022 } 1023 1024 else if(strstr(buffer, " This molecule is "))//rotational data 1025 { 1026 if(strstr(buffer, "asymmetric")) 1027 RotorType = OBRotationData::ASYMMETRIC; 1028 else if(strstr(buffer, "symmetric")) 1029 RotorType = OBRotationData::SYMMETRIC; 1030 else if(strstr(buffer, "linear")) 1031 RotorType = OBRotationData::LINEAR; 1032 else 1033 RotorType = OBRotationData::UNKNOWN; 1034 ifs.getline(buffer,BUFF_SIZE); //symmetry number 1035 tokenize(vs, buffer); 1036 RotSymNum = atoi(vs[3].c_str()); 1037 } 1038 1039 else if(strstr(buffer, "Rotational constant")) 1040 { 1041 tokenize(vs, buffer); 1042 RotConsts.clear(); 1043 for (unsigned int i=3; i<vs.size(); ++i) 1044 RotConsts.push_back(atof(vs[i].c_str())); 1045 } 1046 1047 else if(strstr(buffer, "alpha electrons")) // # of electrons / orbital 1048 { 1049 tokenize(vs, buffer); 1050 if (vs.size() == 6) { 1051 // # alpha electrons # beta electrons 1052 aHOMO = atoi(vs[0].c_str()); 1053 bHOMO = atoi(vs[3].c_str()); 1054 } 1055 } 1056 else if(strstr(buffer, "rbital symmetries")) // orbital symmetries 1057 { 1058 symmetries.clear(); 1059 std::string label; // used as a temporary to remove "(" and ")" from labels 1060 int iii,offset = 0; 1061 bool bDoneSymm; 1062 1063 // Extract both Alpha and Beta symmetries 1064 ifs.getline(buffer, BUFF_SIZE); // skip the current line 1065 for(iii=0; (iii<2); iii++) { 1066 if (strstr(buffer, "electronic state")) 1067 break; // We've gone too far! 1068 while (!ifs.eof() && 1069 (nullptr != strstr(buffer,"Alpha") || 1070 nullptr != strstr(buffer,"Beta"))) { 1071 // skip the Alpha: and Beta: title lines 1072 ifs.getline(buffer, BUFF_SIZE); 1073 } 1074 do { 1075 bDoneSymm = nullptr == strstr(buffer, "("); 1076 if (!bDoneSymm) { 1077 tokenize(vs, buffer); 1078 1079 if (nullptr != strstr(buffer, "Occupied") || nullptr != strstr(buffer, "Virtual")) { 1080 offset = 1; // skip first token 1081 } else { 1082 offset = 0; 1083 } 1084 for (unsigned int i = offset; i < vs.size(); ++i) { 1085 label = vs[i].substr(1, vs[i].length() - 2); 1086 symmetries.push_back(label); 1087 } 1088 ifs.getline(buffer, BUFF_SIZE); // get a new line if we've been reading symmetries 1089 } 1090 // don't read a new line if we're done with symmetries 1091 } while (!ifs.eof() && !bDoneSymm); 1092 } // end alpha/beta section 1093 } 1094 else if (strstr(buffer, "Alpha") && strstr(buffer, ". eigenvalues --")) { 1095 orbitals.clear(); 1096 betaStart = 0; 1097 while (strstr(buffer, ". eigenvalues --")) { 1098 tokenize(vs, buffer); 1099 if (vs.size() < 4) 1100 break; 1101 if (vs[0].find("Beta") !=string::npos && betaStart == 0) // mark where we switch from alpha to beta 1102 betaStart = orbitals.size(); 1103 for (unsigned int i = 4; i < vs.size(); ++i) { 1104 orbitals.push_back(atof(vs[i].c_str())); 1105 } 1106 ifs.getline(buffer, BUFF_SIZE); 1107 } 1108 } 1109 else if(strstr(buffer, " Excited State")) // Force and wavelength data 1110 { 1111 // The above line appears for each state, so just append the info to the vectors 1112 tokenize(vs, buffer); 1113 if (vs.size() >= 9) { 1114 double wavelength = atof(vs[6].c_str()); 1115 double force = atof(vs[8].substr(2).c_str()); // remove the "f=" part 1116 Forces.push_back(force); 1117 Wavelengths.push_back(wavelength); 1118 } 1119 } 1120 else if(strstr(buffer, " Ground to excited state Transition electric dipole moments (Au):")) 1121 // Electronic dipole moments 1122 { 1123 ifs.getline(buffer, BUFF_SIZE); // Headings 1124 ifs.getline(buffer, BUFF_SIZE); // First entry 1125 tokenize(vs, buffer); 1126 while (vs.size() == 5) { 1127 double s = atof(vs[4].c_str()); 1128 EDipole.push_back(s); 1129 ifs.getline(buffer, BUFF_SIZE); 1130 tokenize(vs, buffer); 1131 } 1132 } 1133 else if(strstr(buffer, " state X Y Z R(velocity)")) { 1134 // Rotatory Strengths 1135 ifs.getline(buffer, BUFF_SIZE); // First entry 1136 tokenize(vs, buffer); 1137 while (vs.size() == 5) { 1138 double s = atof(vs[4].c_str()); 1139 RotatoryStrengthsVelocity.push_back(s); 1140 ifs.getline(buffer, BUFF_SIZE); 1141 tokenize(vs, buffer); 1142 } 1143 } 1144 else if(strstr(buffer, " state X Y Z R(length)")) { 1145 // Rotatory Strengths 1146 ifs.getline(buffer, BUFF_SIZE); // First entry 1147 tokenize(vs, buffer); 1148 while (vs.size() == 5) { 1149 double s = atof(vs[4].c_str()); 1150 RotatoryStrengthsLength.push_back(s); 1151 ifs.getline(buffer, BUFF_SIZE); 1152 tokenize(vs, buffer); 1153 } 1154 } 1155 1156 else if (strstr(buffer, "Forces (Hartrees/Bohr)")) 1157 { 1158 ifs.getline(buffer, BUFF_SIZE); // column headers 1159 ifs.getline(buffer, BUFF_SIZE); // ------ 1160 ifs.getline(buffer, BUFF_SIZE); // real data 1161 } 1162 1163 else if (strstr(buffer, "Isotropic = ")) // NMR shifts 1164 { 1165 tokenize(vs, buffer); 1166 if (vs.size() >= 4) 1167 { 1168 atom = mol.GetAtom(atoi(vs[0].c_str())); 1169 OBPairData *nmrShift = new OBPairData(); 1170 nmrShift->SetAttribute("NMR Isotropic Shift"); 1171 1172 string shift = vs[4].c_str(); 1173 nmrShift->SetValue(shift); 1174 1175 atom->SetData(nmrShift); 1176 } 1177 } 1178 else if (strstr(buffer, "SCF Done:") != nullptr) 1179 { 1180 tokenize(vs,buffer); 1181 mol.SetEnergy(atof(vs[4].c_str()) * HARTEE_TO_KCALPERMOL); 1182 confEnergies.push_back(mol.GetEnergy()); 1183 } 1184 /* Temporarily commented out until the handling of energy in OBMol is sorted out 1185 // MP2 energies also use a different syntax 1186 1187 // PM3 energies use a different syntax 1188 else if(strstr(buffer,"E (Thermal)") != nullptr) 1189 { 1190 ifs.getline(buffer,BUFF_SIZE); //Headers 1191 ifs.getline(buffer,BUFF_SIZE); //Total energy; what we want 1192 tokenize(vs,buffer); 1193 mol.SetEnergy(atof(vs[1].c_str())); 1194 confEnergies.push_back(mol.GetEnergy()); 1195 } 1196 */ 1197 else if (strstr(buffer, "Standard basis:") != nullptr) 1198 { 1199 add_unique_pairdata_to_mol(&mol,"basis",buffer,2); 1200 } 1201 else if (strstr(buffer, "Zero-point correction=") != nullptr) 1202 { 1203 tokenize(vs,buffer); 1204 ezpe = atof(vs[2].c_str()); 1205 ezpe_set = true; 1206 } 1207 else if (strstr(buffer, "Thermal correction to Enthalpy=") != nullptr) 1208 { 1209 tokenize(vs,buffer); 1210 Hcorr = atof(vs[4].c_str()); 1211 Hcorr_set = true; 1212 } 1213 else if (strstr(buffer, "Thermal correction to Gibbs Free Energy=") != nullptr) 1214 { 1215 tokenize(vs,buffer); 1216 Gcorr = atof(vs[6].c_str()); 1217 Gcorr_set = true; 1218 } 1219 else if (strstr(buffer, "CV") != nullptr) 1220 { 1221 ifs.getline(buffer,BUFF_SIZE); //Headers 1222 ifs.getline(buffer,BUFF_SIZE); //Total heat capacity 1223 tokenize(vs,buffer); 1224 if (vs.size() == 4) 1225 { 1226 if (vs[0].compare("Total") == 0) 1227 { 1228 CV = atof(vs[2].c_str()); 1229 CV_set = true; 1230 } 1231 } 1232 ifs.getline(buffer,BUFF_SIZE); //Electronic 1233 ifs.getline(buffer,BUFF_SIZE); //Translational 1234 tokenize(vs,buffer); 1235 if ((vs.size() == 4) && (vs[0].compare("Translational") == 0) ) 1236 { 1237 Scomponents.push_back(atof(vs[3].c_str())); 1238 } 1239 ifs.getline(buffer,BUFF_SIZE); //Rotational 1240 tokenize(vs,buffer); 1241 if ((vs.size() == 4) && (vs[0].compare("Rotational") == 0)) 1242 { 1243 Scomponents.push_back(atof(vs[3].c_str())); 1244 } 1245 ifs.getline(buffer,BUFF_SIZE); //Vibrational 1246 tokenize(vs,buffer); 1247 if ((vs.size() == 4) && (vs[0].compare("Vibrational") == 0)) 1248 { 1249 Scomponents.push_back(atof(vs[3].c_str())); 1250 } 1251 } 1252 else if (strstr(buffer,"Temperature=") != nullptr && 1253 strstr(buffer,"Pressure=") != nullptr) 1254 { 1255 tokenize(vs,buffer); 1256 temperature = atof(vs[1].c_str()); 1257 } 1258 else if (strstr(buffer, "(0 K)") != nullptr) 1259 { 1260 /* This must be the last else */ 1261 int i,nsearch; 1262 const char *search[] = { "CBS-QB3 (0 K)", "G2(0 K)", "G3(0 K)", "G4(0 K)", "W1BD (0 K)", "W1U (0 K)" }; 1263 const char *mymeth[] = { "CBS-QB3", "G2", "G3", "G4", "W1BD", "W1U" }; 1264 const int myindex[] = { 3, 2, 2, 2, 3, 3 }; 1265 1266 nsearch = sizeof(search)/sizeof(search[0]); 1267 for(i=0; (i<nsearch); i++) 1268 { 1269 if (strstr(buffer, search[i]) != nullptr) 1270 { 1271 tokenize(vs,buffer); 1272 E0 = atof(vs[myindex[i]].c_str()); 1273 E0_set = 1; 1274 thermo_method = mymeth[i]; 1275 break; 1276 } 1277 } 1278 } 1279 } // end while 1280 1281 if (mol.NumAtoms() == 0) { // e.g., if we're at the end of a file PR#1737209 1282 mol.EndModify(); 1283 return false; 1284 } 1285 1286 mol.EndModify(); 1287 1288 // Check whether we have data to extract heat of formation. 1289 if (ezpe_set && Hcorr_set && Gcorr_set && E0_set && 1290 CV_set && (thermo_method.size() > 0)) 1291 { 1292 extract_thermo(&mol,thermo_method,temperature,ezpe, 1293 Hcorr,Gcorr,E0,CV,RotSymNum,Scomponents); 1294 } 1295 1296 // Attach orbital data, if there is any 1297 if (orbitals.size() > 0) 1298 { 1299 OBOrbitalData *od = new OBOrbitalData; 1300 if (aHOMO == bHOMO) { 1301 od->LoadClosedShellOrbitals(orbitals, symmetries, aHOMO); 1302 } else { 1303 // we have to separate the alpha and beta vectors 1304 std::vector<double> betaOrbitals; 1305 std::vector<std::string> betaSymmetries; 1306 unsigned int initialSize = orbitals.size(); 1307 unsigned int symmSize = symmetries.size(); 1308 if (initialSize != symmSize || betaStart == -1) 1309 { 1310 cerr << "Inconsistency: orbitals have " << initialSize << " elements while symmetries have " << symmSize << endl; 1311 } 1312 else 1313 { 1314 for (unsigned int i = betaStart; i < initialSize; ++i) { 1315 betaOrbitals.push_back(orbitals[i]); 1316 if (symmetries.size() > 0) 1317 betaSymmetries.push_back(symmetries[i]); 1318 } 1319 // ok, now erase the end elements of orbitals and symmetries 1320 for (unsigned int i = betaStart; i < initialSize; ++i) { 1321 orbitals.pop_back(); 1322 if (symmetries.size() > 0) 1323 symmetries.pop_back(); 1324 } 1325 // and load the alphas and betas 1326 od->LoadAlphaOrbitals(orbitals, symmetries, aHOMO); 1327 od->LoadBetaOrbitals(betaOrbitals, betaSymmetries, bHOMO); 1328 } 1329 } 1330 od->SetOrigin(fileformatInput); 1331 mol.SetData(od); 1332 } 1333 1334 //Attach vibrational data, if there is any, to molecule 1335 if(Frequencies.size()>0) 1336 { 1337 OBVibrationData* vd = new OBVibrationData; 1338 vd->SetData(Lx, Frequencies, Intensities); 1339 vd->SetOrigin(fileformatInput); 1340 mol.SetData(vd); 1341 } 1342 //Attach rotational data, if there is any, to molecule 1343 if(RotConsts[0]!=0.0) 1344 { 1345 OBRotationData* rd = new OBRotationData; 1346 rd->SetData(RotorType, RotConsts, RotSymNum); 1347 rd->SetOrigin(fileformatInput); 1348 mol.SetData(rd); 1349 } 1350 // Attach unit cell translation vectors if found 1351 if (numTranslationVectors > 0) { 1352 OBUnitCell* uc = new OBUnitCell; 1353 uc->SetData(translationVectors[0], translationVectors[1], translationVectors[2]); 1354 uc->SetOrigin(fileformatInput); 1355 mol.SetData(uc); 1356 } 1357 //Attach electronic transition data, if there is any, to molecule 1358 if(Forces.size() > 0 && Forces.size() == Wavelengths.size()) 1359 { 1360 OBElectronicTransitionData* etd = new OBElectronicTransitionData; 1361 etd->SetData(Wavelengths, Forces); 1362 if (EDipole.size() == Forces.size()) 1363 etd->SetEDipole(EDipole); 1364 if (RotatoryStrengthsLength.size() == Forces.size()) 1365 etd->SetRotatoryStrengthsLength(RotatoryStrengthsLength); 1366 if (RotatoryStrengthsVelocity.size() == Forces.size()) 1367 etd->SetRotatoryStrengthsVelocity(RotatoryStrengthsVelocity); 1368 etd->SetOrigin(fileformatInput); 1369 mol.SetData(etd); 1370 } 1371 1372 // set some default coordinates 1373 // ConnectTheDots will remove conformers, so we add those later 1374 mol.SetCoordinates(vconf[vconf.size() - 1]); 1375 1376 if (!pConv->IsOption("b",OBConversion::INOPTIONS)) 1377 mol.ConnectTheDots(); 1378 if (!pConv->IsOption("s",OBConversion::INOPTIONS) && !pConv->IsOption("b",OBConversion::INOPTIONS)) 1379 mol.PerceiveBondOrders(); 1380 1381 // Set conformers to all coordinates we adopted 1382 // but remove last geometry -- it's a duplicate 1383 if (vconf.size() > 1) 1384 vconf.pop_back(); 1385 1386 mol.SetConformers(vconf); 1387 mol.SetConformer(mol.NumConformers() - 1); 1388 // Copy the conformer data too 1389 confData->SetDimension(confDimensions); 1390 confData->SetEnergies(confEnergies); 1391 confData->SetForces(confForces); 1392 mol.SetData(confData); 1393 1394 if (hasPartialCharges) { 1395 mol.SetPartialChargesPerceived(); 1396 1397 // Annotate that partial charges come from Mulliken 1398 OBPairData *dp = new OBPairData; 1399 dp->SetAttribute("PartialCharges"); 1400 dp->SetValue(chargeModel); // Mulliken, ESP, etc. 1401 dp->SetOrigin(fileformatInput); 1402 mol.SetData(dp); 1403 } 1404 mol.SetTotalCharge(total_charge); 1405 mol.SetTotalSpinMultiplicity(spin_multiplicity); 1406 1407 mol.SetTitle(title); 1408 return(true); 1409 } 1410 1411 } //namespace OpenBabel 1412