1 /********************************************************************** 2 Copyright (C) 2007 by Maxim Fedorovsky, University of Fribourg (Switzerland). 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 <numeric> 15 #include <typeinfo> 16 #include <functional> 17 #include <cstdlib> 18 #include <algorithm> 19 20 // No diagnoalization yet. Perhaps for 2.2 -GRH 21 // #include <eigen/matrix.h> 22 23 #include <openbabel/mol.h> 24 #include <openbabel/obconversion.h> 25 #include <openbabel/obmolecformat.h> 26 #include <openbabel/mol.h> 27 #include <openbabel/atom.h> 28 #include <openbabel/bond.h> 29 #include <openbabel/obiter.h> 30 #include <openbabel/elements.h> 31 #include <openbabel/generic.h> 32 33 #include <openbabel/math/matrix3x3.h> 34 35 #define BOHR2ANGSTROM 0.5291772083 36 #define HARTREE2INVCM 219474.631371 37 #define AMU2AU 1822.88848 38 #define REDDIPSTR2INT 974.864 39 40 using namespace std; 41 42 namespace OpenBabel 43 { 44 //! \brief Class for parsing Gaussian formatted checkpoint files 45 class FCHKFormat : public OBMoleculeFormat 46 { 47 public : 48 /* Register this format type ID */ FCHKFormat()49 FCHKFormat() 50 { 51 OBConversion::RegisterFormat("fchk", this, 52 "chemical/x-gaussian-checkpoint"); 53 OBConversion::RegisterFormat("fch", this, 54 "chemical/x-gaussian-checkpoint"); 55 OBConversion::RegisterFormat("fck", this, 56 "chemical/x-gaussian-checkpoint"); 57 } 58 Description()59 virtual const char * Description() 60 { 61 return "Gaussian formatted checkpoint file format\n" 62 "A formatted text file containing the results of a Gaussian calculation\n" 63 "Currently supports reading molecular geometries from fchk files. More to come.\n\n" 64 65 "Read options e.g. -as\n" 66 " s Single bonds only\n" 67 " b No bond perception\n\n"; 68 // Vibrational analysis not yet supported in OB-2.1. 69 // v Do not perform the vibrational analysis\n\n"; 70 }; 71 GetMIMEType()72 virtual const char * GetMIMEType() 73 { 74 return "chemical/x-gaussian-checkpoint"; 75 }; 76 Flags()77 virtual unsigned int Flags() 78 { 79 return NOTWRITABLE; 80 }; 81 82 virtual bool ReadMolecule(OBBase *, OBConversion *); 83 84 private : 85 bool static read_int(const char * const, int * const); 86 template<class T> static bool read_numbers(const char * const, 87 vector<T> &, 88 const unsigned int width = 0); 89 template <class T> static bool read_section(const char * const, 90 vector<T> &, 91 const unsigned int, 92 bool * const, 93 const char * const, 94 const unsigned int, 95 const unsigned int width = 0); 96 bool static validate_section(const char * const, 97 const int, 98 const char * const, 99 const unsigned int); 100 bool static validate_number(const int, 101 const char * const, 102 const unsigned int); 103 void static construct_mol(OBMol * const, 104 OBConversion * const, 105 const unsigned int, 106 const vector<int> &, 107 const vector<double> &, 108 const int, 109 const vector<int> &, 110 const vector<int> &); 111 // void static vibana(OBMol * const, 112 // const vector<double> &, 113 // const vector<double> &); 114 // bool static generate_tr_rot(OBMol * const, 115 // double * const, 116 // unsigned int * const); 117 // void static rotmatrix_axis(const double * const, const double, 118 // matrix3x3 &); 119 // double static cosine(const double * const, const double * const); 120 // void static normalize_set(double * const, 121 // const unsigned int, const unsigned int); 122 // void static orthogonalize_set(double * const, 123 // const unsigned int, const unsigned int, 124 // const bool); 125 // void static rests(const unsigned int, unsigned int * const, 126 // unsigned int * const, const unsigned int); 127 128 // void static retrieve_hessian_indices(const unsigned int, 129 // unsigned int * const, 130 // unsigned int * const, 131 // unsigned int * const, 132 // unsigned int * const); 133 }; 134 135 /* Make an instance of the format class */ 136 FCHKFormat theFCHKFormat; 137 138 139 /*! 140 **\brief Extract the data. 141 ** 142 **The data are : 143 ** comment 144 ** charge 145 ** multiplicity 146 ** geometry 147 ** connection table 148 */ ReadMolecule(OBBase * pOb,OBConversion * pConv)149 bool FCHKFormat::ReadMolecule(OBBase * pOb, OBConversion * pConv) 150 { 151 OBMol * const pmol = dynamic_cast<OBMol*>(pOb); 152 153 if (!pmol) 154 return false; 155 156 /* input stream */ 157 istream * const pifs = pConv->GetInStream(); 158 159 /* for logging errors */ 160 stringstream error_msg; 161 162 /* help variables */ 163 char buff[BUFF_SIZE]; 164 unsigned int lineno = 0; 165 int intvar; 166 bool finished; 167 bool atomnos_found = false; 168 bool coords_found = false; 169 bool nbond_found = false; 170 bool ibond_found = false; 171 bool hessian_found = false; 172 bool dipder_found = false; 173 bool alphaorb_found = false; 174 bool betaorb_found = false; 175 176 /* variables for storing the read data */ 177 int Natoms = -1, MxBond = -1; 178 int numAOrb = -1, numBOrb = -1; 179 int numAElec = -1, numBElec = -1; 180 vector<int> atomnos, NBond, IBond; 181 vector<double> coords, hessian, dipder, alphaorb, betaorb; 182 183 /*** start reading ***/ 184 pmol->BeginModify(); 185 186 /* first line : comment */ 187 if (!pifs->getline(buff, BUFF_SIZE)) 188 { 189 obErrorLog.ThrowError("FCHKFormat::ReadMolecule()", 190 "Failed to read the comment.", 191 obError); 192 return false; 193 } 194 195 ++lineno; 196 pmol->SetTitle(buff); 197 198 /* extract all the information here */ 199 while (pifs->getline(buff, BUFF_SIZE)) 200 { 201 ++lineno; 202 203 /* Number of atoms */ 204 if (buff == strstr(buff, "Number of atoms")) 205 { 206 if (!FCHKFormat::read_int(buff, &Natoms)) 207 { 208 error_msg << "Could not read the number of atoms from line #" 209 << lineno << "."; 210 211 obErrorLog.ThrowError("FCHKFormat::ReadMolecule()", 212 error_msg.str(), 213 obError); 214 return false; 215 } 216 if (0 >= Natoms) 217 { 218 error_msg << "Invalid number of atoms : " << Natoms << "."; 219 obErrorLog.ThrowError("FCHKFormat::ReadMolecule()", 220 error_msg.str(), 221 obError); 222 return false; 223 } 224 continue; 225 } 226 227 /* Dipole Moment */ 228 if (buff == strstr(buff, "Dipole Moment")) 229 { 230 if (!pifs->getline(buff, BUFF_SIZE)) { 231 error_msg << "Could not read the dipole moment after line #" 232 << lineno << "."; 233 obErrorLog.ThrowError("FCHKFormat::ReadMolecule()", 234 error_msg.str(), 235 obError); 236 return false; 237 } 238 ++lineno; 239 vector<double> moments; 240 if (!FCHKFormat::read_numbers(buff, moments) || moments.size() != 3) { 241 error_msg << "Could not read the dipole moment from line #" 242 << lineno << "."; 243 obErrorLog.ThrowError("FCHKFormat::ReadMolecule()", 244 error_msg.str(), 245 obError); 246 return false; 247 } 248 249 OBVectorData *dipoleMoment = new OBVectorData; 250 dipoleMoment->SetAttribute("Dipole Moment"); 251 dipoleMoment->SetData(moments[0], moments[1], moments[2]); 252 dipoleMoment->SetOrigin(fileformatInput); 253 pmol->SetData(dipoleMoment); 254 continue; 255 } 256 257 /* Charge */ 258 if (buff == strstr(buff, "Charge")) 259 { 260 if (!FCHKFormat::read_int(buff, &intvar)) 261 { 262 error_msg << "Could not read the charge from line #" 263 << lineno << "."; 264 obErrorLog.ThrowError("FCHKFormat::ReadMolecule()", 265 error_msg.str(), 266 obError); 267 return false; 268 } 269 270 pmol->SetTotalCharge(intvar); 271 continue; 272 } 273 274 /* Multiplicity */ 275 if (buff == strstr(buff, "Multiplicity")) 276 { 277 if (!FCHKFormat::read_int(buff, &intvar)) 278 { 279 error_msg << "Could not read the multiplicity from line #" 280 << lineno << "."; 281 obErrorLog.ThrowError("FCHKFormat::ReadMolecule()", 282 error_msg.str(), 283 obError); 284 return false; 285 } 286 287 pmol->SetTotalSpinMultiplicity(intvar); 288 continue; 289 } 290 291 /* start of the "Atomic numbers" section */ 292 if (buff == strstr(buff, "Atomic numbers")) 293 { 294 if (!FCHKFormat::validate_number(Natoms, "number of atoms", lineno)) 295 return false; 296 297 if (!FCHKFormat::validate_section(buff, 298 Natoms, 299 "number of atomic numbers", 300 lineno)) 301 return false; 302 303 atomnos_found = true; 304 continue; 305 } 306 307 /* start of the "Current cartesian coordinates" section */ 308 if (buff == strstr(buff, "Current cartesian coordinates")) 309 { 310 if (!FCHKFormat::validate_number(Natoms, "number of atoms", lineno)) 311 return false; 312 313 if (!FCHKFormat::validate_section(buff, 314 3 * Natoms, 315 "number of coordinates", 316 lineno)) 317 return false; 318 319 coords_found = true; 320 continue; 321 } 322 323 /* MxBond - largest number of bonds to an atom */ 324 if (buff == strstr(buff, "MxBond")) 325 { 326 if (!FCHKFormat::read_int(buff, &MxBond)) 327 { 328 error_msg << "Could not read the number of bonds to each atom" 329 << " (MxBond) from line #" << lineno << "."; 330 331 obErrorLog.ThrowError("FCHKFormat::ReadMolecule()", 332 error_msg.str(), 333 obError); 334 return false; 335 } 336 if (0 >= MxBond) 337 { 338 error_msg << "Invalid number of bonds to each atom : " 339 << MxBond << "."; 340 obErrorLog.ThrowError("FCHKFormat::ReadMolecule()", 341 error_msg.str(), 342 obError); 343 return false; 344 } 345 continue; 346 } 347 348 /* start of the "NBond" section */ 349 if (buff == strstr(buff, "NBond")) 350 { 351 if (!FCHKFormat::validate_number(Natoms, "number of atoms", lineno)) 352 return false; 353 354 if (!FCHKFormat::validate_number(MxBond, 355 "number of bonds to each atom", 356 lineno)) 357 return false; 358 359 if (!FCHKFormat::validate_section(buff, 360 Natoms, 361 "number of bonds to each atom", 362 lineno)) 363 return false; 364 365 nbond_found = true; 366 continue; 367 } 368 369 /* start of the "IBond" section */ 370 if (buff == strstr(buff, "IBond")) 371 { 372 if (!FCHKFormat::validate_number(Natoms, "number of atoms", lineno)) 373 return false; 374 375 if (!FCHKFormat::validate_number(MxBond, 376 "number of bonds to each atom", 377 lineno)) 378 return false; 379 380 if (!FCHKFormat::validate_section(buff, 381 MxBond * Natoms, 382 "number of bonds", 383 lineno)) 384 return false; 385 386 ibond_found = true; 387 continue; 388 } 389 390 /* start of the "Cartesian Force Constants" section */ 391 if (buff == strstr(buff, "Cartesian Force Constants")) 392 { 393 if (!FCHKFormat::validate_number(Natoms, "number of atoms", lineno)) 394 return false; 395 396 if (!FCHKFormat::validate_section(buff, 397 (3 * Natoms) * (3 * Natoms + 1) / 2, 398 "number of force constants", 399 lineno)) 400 return false; 401 402 hessian_found = true; 403 continue; 404 } 405 406 /* start of the "Dipole Derivatives" section */ 407 if (buff == strstr(buff, "Dipole Derivatives")) 408 { 409 if (!FCHKFormat::validate_number(Natoms, "number of atoms", lineno)) 410 return false; 411 412 if (!FCHKFormat::validate_section(buff, 413 9 * Natoms, 414 "number of dipole derivatives", 415 lineno)) 416 return false; 417 418 dipder_found = true; 419 continue; 420 } 421 422 if (buff == strstr(buff, "Number of alpha electrons")) 423 { 424 if (!FCHKFormat::read_int(buff, &numAElec)) 425 { 426 error_msg << "Could not read the number of alpha electrons" 427 << " from line #" << lineno << "."; 428 429 obErrorLog.ThrowError("FCHKFormat::ReadMolecule()", 430 error_msg.str(), 431 obError); 432 return false; 433 } 434 if (0 >= numAElec) 435 { 436 error_msg << "Invalid number of alpha electrons: " 437 << MxBond << "."; 438 obErrorLog.ThrowError("FCHKFormat::ReadMolecule()", 439 error_msg.str(), 440 obError); 441 return false; 442 } 443 continue; 444 } 445 446 if (buff == strstr(buff, "Number of beta electrons")) 447 { 448 if (!FCHKFormat::read_int(buff, &numBElec)) 449 { 450 error_msg << "Could not read the number of beta electrons" 451 << " from line #" << lineno << "."; 452 453 obErrorLog.ThrowError("FCHKFormat::ReadMolecule()", 454 error_msg.str(), 455 obError); 456 return false; 457 } 458 if (0 >= numBElec) 459 { 460 error_msg << "Invalid number of beta electrons: " 461 << MxBond << "."; 462 obErrorLog.ThrowError("FCHKFormat::ReadMolecule()", 463 error_msg.str(), 464 obError); 465 return false; 466 } 467 continue; 468 } 469 470 /* start of the alpha orbitals section */ 471 if (buff == strstr(buff, "Alpha Orbital Energies")) 472 { 473 if (!FCHKFormat::read_int(buff, &numAOrb)) 474 { 475 error_msg << "Could not read the number of alpha orbitals" 476 << " from line #" << lineno << "."; 477 478 obErrorLog.ThrowError("FCHKFormat::ReadMolecule()", 479 error_msg.str(), 480 obError); 481 return false; 482 } 483 if (0 >= numAOrb) 484 { 485 error_msg << "Invalid number of alpha orbitals: " 486 << MxBond << "."; 487 obErrorLog.ThrowError("FCHKFormat::ReadMolecule()", 488 error_msg.str(), 489 obError); 490 return false; 491 } 492 493 alphaorb_found = true; 494 continue; 495 } 496 497 /* start of the beta orbitals section */ 498 if (buff == strstr(buff, "Beta Orbital Energies")) 499 { 500 if (!FCHKFormat::read_int(buff, &numBOrb)) 501 { 502 error_msg << "Could not read the number of beta orbitals" 503 << " from line #" << lineno << "."; 504 505 obErrorLog.ThrowError("FCHKFormat::ReadMolecule()", 506 error_msg.str(), 507 obError); 508 return false; 509 } 510 if (0 >= numBOrb) 511 { 512 error_msg << "Invalid number of beta orbitals: " 513 << MxBond << "."; 514 obErrorLog.ThrowError("FCHKFormat::ReadMolecule()", 515 error_msg.str(), 516 obError); 517 return false; 518 } 519 520 betaorb_found = true; 521 continue; 522 } 523 524 /* reading the atomic numbers */ 525 if (atomnos_found && Natoms > (int)atomnos.size()) 526 { 527 if (!FCHKFormat::read_section(buff, atomnos, Natoms, &finished, 528 "atomic numbers", lineno)) 529 return false; 530 531 atomnos_found = !finished; 532 continue; 533 } 534 535 /* reading the coordinates */ 536 if (coords_found && 3 * Natoms > (int)coords.size()) 537 { 538 if (!FCHKFormat::read_section(buff, coords, 3 * Natoms, &finished, 539 "coordinates", lineno, 16)) 540 return false; 541 542 coords_found = !finished; 543 continue; 544 } 545 546 /* reading the numbers of bonds to each atom */ 547 if (nbond_found && Natoms > (int)NBond.size()) 548 { 549 if (!FCHKFormat::read_section(buff, NBond, Natoms, &finished, 550 "number of bonds to each atom", lineno)) 551 return false; 552 553 nbond_found = !finished; 554 continue; 555 } 556 557 /* reading the list of atoms bound to each atom */ 558 if (ibond_found && MxBond * Natoms > (int)IBond.size()) 559 { 560 if (!FCHKFormat::read_section(buff, IBond, MxBond * Natoms, &finished, 561 "atom bonds", lineno)) 562 return false; 563 564 ibond_found = !finished; 565 continue; 566 } 567 568 /* reading the hessian */ 569 if (hessian_found && 570 (3 * Natoms) * (3 * Natoms + 1) / 2 > (int)hessian.size()) 571 { 572 if (!FCHKFormat::read_section(buff, hessian, 573 (3 * Natoms) * (3 * Natoms + 1) / 2, 574 &finished, 575 "hessian", lineno)) 576 return false; 577 578 hessian_found = !finished; 579 continue; 580 } 581 582 /* reading the dipole derivatives */ 583 if (dipder_found && 9 * Natoms > (int)dipder.size()) 584 { 585 if (!FCHKFormat::read_section(buff, dipder, 586 9 * Natoms, 587 &finished, 588 "dipole derivatives", lineno)) 589 return false; 590 591 dipder_found = !finished; 592 continue; 593 } 594 595 /* reading the alpha orbital energies */ 596 if (alphaorb_found && (numAOrb > alphaorb.size())) 597 { 598 if (!FCHKFormat::read_section(buff, alphaorb, 599 numAOrb, 600 &finished, 601 "alpha orbital energies", lineno)) 602 return false; 603 604 alphaorb_found = !finished; // if we've read once, don't re-read 605 continue; 606 } 607 608 /* reading the beta orbital energies */ 609 if (betaorb_found && (numBOrb > betaorb.size())) 610 { 611 if (!FCHKFormat::read_section(buff, betaorb, 612 numBOrb, 613 &finished, 614 "beta orbital energies", lineno)) 615 return false; 616 617 betaorb_found = !finished; 618 continue; 619 } 620 621 } /* while */ 622 623 /*** validating ***/ 624 if (0 >= Natoms) 625 { 626 error_msg << "Number of atoms could not be read."; 627 obErrorLog.ThrowError("FCHKFormat::ReadMolecule()", 628 error_msg.str(), 629 obError); 630 return false; 631 } 632 633 if (atomnos.empty()) 634 { 635 error_msg << "\"Atomic numbers\" section was not found."; 636 obErrorLog.ThrowError("FCHKFormat::ReadMolecule()", 637 error_msg.str(), 638 obError); 639 return false; 640 } 641 642 if (coords.empty()) 643 { 644 error_msg << "\"Current cartesian coordinates\" section was not found."; 645 obErrorLog.ThrowError("FCHKFormat::ReadMolecule()", 646 error_msg.str(), 647 obError); 648 return false; 649 } 650 651 /* connectivity : if MxBond is set, all the sections must be supplied */ 652 if (-1 != MxBond) 653 { 654 if (NBond.empty() || IBond.empty()) 655 { 656 error_msg << "If MxBond is set, then the \"NBond\" and \"IBond\"" 657 << " sections *must* be supplied."; 658 obErrorLog.ThrowError("FCHKFormat::ReadMolecule()", 659 error_msg.str(), 660 obError); 661 return false; 662 } 663 664 /* not nbonds > MxBond or <= 0 665 no atom numbers < 0 or > Natoms */ 666 if (NBond.end() != find_if(NBond.begin(), 667 NBond.end(), 668 bind2nd(less_equal<int>(), 0)) || 669 NBond.end() != find_if(NBond.begin(), 670 NBond.end(), 671 bind2nd(greater<int>(), MxBond)) || 672 IBond.end() != find_if(IBond.begin(), 673 IBond.end(), 674 bind2nd(less<int>(), 0)) || 675 IBond.end() != find_if(IBond.begin(), 676 IBond.end(), 677 bind2nd(greater<int>(), Natoms))) 678 { 679 error_msg << "Invalid connectivity : check the \"NBond\" and/or" 680 << " \"IBond\" section(s)."; 681 obErrorLog.ThrowError("FCHKFormat::ReadMolecule()", 682 error_msg.str(), 683 obWarning); 684 MxBond = -1; 685 } 686 } 687 688 /*** finalizing ***/ 689 FCHKFormat::construct_mol(pmol, pConv, Natoms, atomnos, coords, 690 MxBond, NBond, IBond); 691 692 // if (!hessian.empty() && !pConv->IsOption("v", OBConversion::INOPTIONS)) 693 // { 694 // FCHKFormat::vibana(pmol, hessian, dipder); 695 // } 696 697 /*** finished ***/ 698 pmol->EndModify(); 699 700 if (numAOrb > 0 && alphaorb.size() == numAOrb) { 701 OBOrbitalData *od = new OBOrbitalData; // create new store 702 vector<string> symmetries; 703 704 if (numBOrb > 0 && betaorb.size() == numBOrb) { // open shell calculation 705 od->LoadAlphaOrbitals(alphaorb, symmetries, numAElec); 706 od->LoadBetaOrbitals(betaorb, symmetries, numBElec); 707 } else { 708 od->LoadClosedShellOrbitals(alphaorb, symmetries, numAElec); 709 } 710 od->SetOrigin(fileformatInput); 711 pmol->SetData(od); 712 } 713 714 return true; 715 } 716 717 /*! 718 **\brief Convert the last token in a string to integer 719 ** 720 **Return false if the conversion failed. 721 **The read integer is stored in the num argument. 722 */ read_int(const char * const line,int * const num)723 bool FCHKFormat::read_int(const char * const line, int * const num) 724 { 725 vector<string> vs; 726 char * endptr; 727 728 tokenize(vs, line); 729 *num = strtol(vs.back().c_str(), &endptr, 10); 730 731 return endptr != vs.back().c_str(); 732 } 733 734 /*! 735 **\brief Read integers or doubles into a vector from a string 736 **\param line string 737 **\param v vector to read in 738 */ 739 template<class T> read_numbers(const char * const line,vector<T> & v,const unsigned int width)740 bool FCHKFormat::read_numbers(const char * const line, vector<T> & v, const unsigned int width) 741 { 742 if (width == 0) { 743 vector<string> vs; 744 tokenize(vs, line); 745 746 if (0 < vs.size()) 747 { 748 char * endptr; 749 T num; 750 751 for (vector<string>::const_iterator it=vs.begin(); vs.end() != it; ++it) 752 { 753 if (typeid(double) == typeid(T)) 754 num = static_cast<T>(strtod((*it).c_str(), &endptr)); 755 else 756 num = static_cast<T>(strtol((*it).c_str(), &endptr, 10)); 757 758 /* quit if the value cannot be read */ 759 if (endptr == (*it).c_str()) 760 return false; 761 762 v.push_back(num); 763 } 764 } 765 } else { 766 // fixed width records (e.g., coordinates = 16 characters) 767 string lineCopy(line), substring; 768 T num; 769 char * endptr; 770 int maxColumns = 80 / width; 771 772 for (int column = 0; column < maxColumns; ++column) { 773 substring = lineCopy.substr(column * width, width); 774 775 if (typeid(double) == typeid(T)) 776 num = static_cast<T>(strtod(substring.c_str(), &endptr)); 777 else 778 num = static_cast<T>(strtol(substring.c_str(), &endptr, 10)); 779 780 /* quit if the value cannot be read */ 781 if (endptr == substring.c_str()) 782 break; 783 784 v.push_back(num); 785 } 786 } 787 788 return true; 789 } 790 791 /*! 792 **\brief Validate a section 793 ** 794 **Try to convert the last token of a string to integer and compare it to a 795 **given number. If they are not equal report an error and return false. 796 **Otherwise return true. 797 */ validate_section(const char * const line,const int nreq,const char * const desc,const unsigned int lineno)798 bool FCHKFormat::validate_section(const char * const line, 799 const int nreq, 800 const char * const desc, 801 const unsigned int lineno) 802 { 803 int intvar; 804 stringstream error_msg; 805 806 if (!FCHKFormat::read_int(line, &intvar)) 807 { 808 error_msg << "Could not read the " << desc 809 << " from line #" << lineno << "."; 810 obErrorLog.ThrowError("FCHKFormat::ReadMolecule()", 811 error_msg.str(), 812 obError); 813 return false; 814 } 815 if (intvar != nreq) 816 { 817 error_msg << desc << " must be exactly " << nreq 818 << ", found " << intvar << "."; 819 obErrorLog.ThrowError("FCHKFormat::ReadMolecule()", 820 error_msg.str(), 821 obError); 822 return false; 823 } 824 825 return true; 826 } 827 828 /*! 829 **\brief Validate a number 830 ** 831 **Compare it with a given value and report an error unless they are equal. 832 */ validate_number(const int num,const char * const desc,const unsigned int lineno)833 bool FCHKFormat::validate_number(const int num, 834 const char * const desc, 835 const unsigned int lineno) 836 { 837 stringstream error_msg; 838 839 if (-1 == num) 840 { 841 error_msg << desc << " must be already read before line #" 842 << lineno << "."; 843 obErrorLog.ThrowError("FCHKFormat::ReadMolecule()", 844 error_msg.str(), 845 obError); 846 return false; 847 } 848 return true; 849 } 850 851 /*! 852 **\brief Read numbers in a section 853 **\param line line 854 **\param v vector to read in 855 **\param nreq maximum number of elements in v 856 **\param finished whether the reading is finished 857 **\param desc name of the section 858 **\param lineno line number 859 */ 860 template <class T> read_section(const char * const line,vector<T> & v,const unsigned int nreq,bool * const finished,const char * const desc,const unsigned int lineno,const unsigned int width)861 bool FCHKFormat::read_section(const char * const line, 862 vector<T> & v, 863 const unsigned int nreq, 864 bool * const finished, 865 const char * const desc, 866 const unsigned int lineno, 867 const unsigned int width) 868 { 869 stringstream error_msg; 870 871 *finished = false; 872 873 if(!FCHKFormat::read_numbers(line, v, width)) 874 { 875 error_msg << "Expecting " << desc << " in line #" << lineno << "."; 876 obErrorLog.ThrowError("FCHKFormat::ReadMolecule()", 877 error_msg.str(), 878 obError); 879 return false; 880 } 881 if (nreq <= v.size()) 882 { 883 *finished = true; 884 if (nreq < v.size()) 885 { 886 error_msg << "Ignoring the superfluous " << desc 887 << "in line #" << lineno << "."; 888 obErrorLog.ThrowError("FCHKFormat::ReadMolecule()", 889 error_msg.str(), 890 obWarning); 891 } 892 } 893 894 return true; 895 } 896 897 /*! 898 **\brief Construct the molecule 899 **\param pmol molecule 900 **\param Natoms number of atoms 901 **\param atomnos atomic numbers 902 **\param coords coordinates 903 **\param MxBond largest number of bonds to an atom (usually 4) 904 **\param NBond number of bonds to each atom 905 **\param IBond list of atoms bounded to each atom 906 */ construct_mol(OBMol * const pmol,OBConversion * const pConv,const unsigned int Natoms,const vector<int> & atomnos,const vector<double> & coords,const int MxBond,const vector<int> & NBond,const vector<int> & IBond)907 void FCHKFormat::construct_mol(OBMol * const pmol, 908 OBConversion * const pConv, 909 const unsigned int Natoms, 910 const vector<int> & atomnos, 911 const vector<double> & coords, 912 const int MxBond, 913 const vector<int> & NBond, 914 const vector<int> & IBond) 915 { 916 pmol->ReserveAtoms(Natoms); 917 918 OBAtom * atom; 919 for (unsigned int a = 0; Natoms > a; ++a) 920 { 921 atom = pmol->NewAtom(); 922 923 atom->SetAtomicNum(atomnos[a]); 924 atom->SetVector(coords[0 + 3 * a] * BOHR2ANGSTROM, 925 coords[1 + 3 * a] * BOHR2ANGSTROM, 926 coords[2 + 3 * a] * BOHR2ANGSTROM); 927 } 928 929 /* unless suppressed */ 930 if (!pConv->IsOption("b", OBConversion::INOPTIONS)) 931 { 932 /* if the connectivity is provided */ 933 if (-1 != MxBond) 934 { 935 /* atoms */ 936 for (unsigned int a = 0; Natoms > a; ++a) 937 { 938 /* bonds for an atom */ 939 for (unsigned int i = 0; (unsigned int)NBond[a] > i; ++i) 940 { 941 pmol->AddBond(1 + a, IBond[MxBond*a + i], 1); 942 } 943 } 944 } 945 else 946 { 947 pmol->ConnectTheDots(); 948 } 949 } 950 951 if (!pConv->IsOption("s", OBConversion::INOPTIONS) && 952 !pConv->IsOption("b", OBConversion::INOPTIONS)) 953 pmol->PerceiveBondOrders(); 954 } 955 956 } /* namespace OpenBabel */ 957