1 /********************************************************************** 2 Copyright (C) 2004 by Chris Morley for template 3 Copyright (C) 2009 by David C. Lonie for VASP 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 15 #include <openbabel/babelconfig.h> 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 24 25 #include <limits> 26 #include <locale> // For isalpha(int) 27 #include <functional> 28 #include <iostream> 29 #include <algorithm> 30 #include <cstdlib> 31 32 #define EV_TO_KCAL_PER_MOL 23.060538 33 34 using namespace std; 35 namespace OpenBabel { 36 class VASPFormat : public OBMoleculeFormat 37 { 38 protected: 39 class compare_sort_items 40 { 41 std::vector<int> csm; 42 bool num_sort; 43 public: compare_sort_items(const std::vector<int> & _custom_sort_nums,bool _num_sort)44 compare_sort_items(const std::vector<int> &_custom_sort_nums, bool _num_sort): 45 csm(_custom_sort_nums), num_sort(_num_sort) {}; operator ()(const OBAtom * a,const OBAtom * b)46 bool operator()(const OBAtom *a, const OBAtom *b) 47 { 48 int a_num = a->GetAtomicNum(); 49 int b_num = b->GetAtomicNum(); 50 int dist = std::distance(std::find(csm.begin(), csm.end(), b_num), 51 std::find(csm.begin(), csm.end(), a_num)); 52 53 if ( dist != 0) 54 return dist < 0; 55 56 if( (num_sort) && ( a_num - b_num != 0 ) ) 57 return a_num < b_num; 58 59 return false; 60 } 61 }; 62 public: 63 VASPFormat()64 VASPFormat() 65 { 66 // This will actually read the CONTCAR file: 67 OBConversion::RegisterFormat("CONTCAR",this); 68 OBConversion::RegisterFormat("POSCAR",this); 69 OBConversion::RegisterFormat("VASP",this); 70 OBConversion::RegisterOptionParam("s", this, 0, OBConversion::INOPTIONS); 71 OBConversion::RegisterOptionParam("b", this, 0, OBConversion::INOPTIONS); 72 OBConversion::RegisterOptionParam("w", this, 0, OBConversion::OUTOPTIONS); 73 OBConversion::RegisterOptionParam("z", this, 0, OBConversion::OUTOPTIONS); 74 OBConversion::RegisterOptionParam("4", this, 0, OBConversion::OUTOPTIONS); 75 } 76 Description()77 virtual const char* Description() 78 { 79 return 80 "VASP format\n" 81 "Reads in data from POSCAR and CONTCAR to obtain information from " 82 "VASP calculations.\n\n" 83 84 "Due to limitations in Open Babel's file handling, reading in VASP\n" 85 "files can be a bit tricky; the client that is using Open Babel must\n" 86 "use OBConversion::ReadFile() to begin the conversion. This change is\n" 87 "usually trivial. Also, the complete path to the CONTCAR/POSCAR file\n" 88 "must be provided, otherwise the other files needed will not be\n" 89 "found.\n\n" 90 91 "Both VASP 4.x and 5.x POSCAR formats are supported.\n\n" 92 93 "By default, atoms are written out in the order they are present in the input\n" 94 "molecule. To sort by atomic number specify ``-xw``. To specify the sort\n" 95 "order, use the ``-xz`` option.\n\n" 96 97 "Read Options e.g. -as\n" 98 " s Output single bonds only\n" 99 " b Disable bonding entirely\n\n" 100 101 "Write Options e.g. -x4\n" 102 " w Sort atoms by atomic number\n" 103 " z <list of atoms> Specify the order to write out atoms\n" 104 " 'atom1 atom2 ...': atom1 first, atom2 second, etc. The remaining\n" 105 " atoms are written in the default order or (if ``-xw`` is specified)\n" 106 " in order of atomic number.\n" 107 " 4 Write a POSCAR using the VASP 4.x specification.\n" 108 " The default is to use the VASP 5.x specification.\n\n" 109 ; 110 111 }; 112 SpecificationURL()113 virtual const char* SpecificationURL(){return "http://cms.mpi.univie.ac.at/vasp/vasp/vasp.html";}; 114 115 /* Flags() can return be any of the following combined by | 116 or be omitted if none apply 117 NOTREADABLE READONEONLY NOTWRITABLE WRITEONEONLY DEFAULTFORMAT 118 READBINARY WRITEBINARY READXML ZEROATOMSOK */ Flags()119 virtual unsigned int Flags() 120 { 121 return READONEONLY; 122 }; 123 SkipObjects(int n,OBConversion * pConv)124 virtual int SkipObjects(int n, OBConversion* pConv) 125 { 126 return 0; 127 }; 128 129 //////////////////////////////////////////////////// 130 /// Declarations for the "API" interface functions. Definitions are below 131 virtual bool ReadMolecule(OBBase* pOb, OBConversion* pConv); 132 virtual bool WriteMolecule(OBBase* pOb, OBConversion* pConv); 133 134 private: 135 /* Add declarations for any local function or member variables used. 136 Generally only a single instance of a format class is used. Keep this in 137 mind if you employ member variables. */ 138 }; 139 //////////////////////////////////////////////////// 140 141 //Make an instance of the format class 142 VASPFormat theVASPFormat; 143 144 ///////////////////////////////////////////////////////////////// 145 ReadMolecule(OBBase * pOb,OBConversion * pConv)146 bool VASPFormat::ReadMolecule(OBBase* pOb, OBConversion* pConv) 147 { 148 OBMol* pmol = pOb->CastAndClear<OBMol>(); 149 if (pmol == nullptr) 150 return false; 151 152 // Move stream to EOF, some apps check ifs position to check for multimolecule files. 153 // VASP does not support this, and this parser makes its own streams. 154 istream &ifs = *pConv->GetInStream(); 155 ifs.seekg (0, ios::end); 156 157 char buffer[BUFF_SIZE], tag[BUFF_SIZE]; 158 double x,y,z, scale; 159 unsigned int totalAtoms=0, atomIndex=0, atomCount=0; 160 OBAtom *atom; 161 bool cartesian; 162 string str, path; 163 vector<string> vs; 164 vector<unsigned int> numAtoms, atomTypes; 165 bool selective; // is selective dynamics used? 166 string key, value; // store the info about constraints 167 OBPairData *cp; // in this PairData 168 bool hasEnthalpy=false; 169 bool hasVibrations=false; 170 bool needSymbolsInGeometryFile = false; 171 double enthalpy_eV, pv_eV; 172 vector<vector <vector3> > Lx; 173 vector<double> Frequencies; 174 vector<matrix3x3> dipGrad; 175 176 // Get path of CONTCAR/POSCAR: 177 // ifs_path.getline(buffer,BUFF_SIZE); 178 // path = buffer; 179 path = pConv->GetInFilename(); 180 if (path.empty()) return false; // Should be using ReadFile, not Read! 181 size_t found; 182 found = path.rfind("/"); 183 path = path.substr(0, found); 184 if (found == string::npos) path = "./"; // No "/" in path? 185 186 // Open files 187 string potcar_filename = path + "/POTCAR"; 188 string outcar_filename = path + "/OUTCAR"; 189 string doscar_filename = path + "/DOSCAR"; 190 string contcar_filename = pConv->GetInFilename(); // POSCAR _OR_ CONTCAR 191 ifstream ifs_pot (potcar_filename.c_str()); 192 ifstream ifs_out (outcar_filename.c_str()); 193 ifstream ifs_dos (doscar_filename.c_str()); 194 ifstream ifs_cont (contcar_filename.c_str()); 195 if (!ifs_pot) { 196 needSymbolsInGeometryFile = true; 197 } 198 if (!ifs_cont) { 199 return false; // No geometry file? 200 } 201 202 pmol->BeginModify(); 203 204 // Start working on CONTCAR: 205 ifs_cont.getline(buffer,BUFF_SIZE); // Comment line 206 pmol->SetTitle(buffer); 207 ifs_cont.getline(buffer,BUFF_SIZE); // Scale 208 scale = atof(buffer); 209 210 ifs_cont.getline(buffer,BUFF_SIZE); // X_Vec vector 211 tokenize(vs, buffer); 212 x = atof(vs.at(0).c_str()) * scale; 213 y = atof(vs.at(1).c_str()) * scale; 214 z = atof(vs.at(2).c_str()) * scale; 215 vector3 x_vec (x,y,z); 216 217 ifs_cont.getline(buffer,BUFF_SIZE); // Y_Vec vector 218 tokenize(vs, buffer); 219 x = atof(vs.at(0).c_str()) * scale; 220 y = atof(vs.at(1).c_str()) * scale; 221 z = atof(vs.at(2).c_str()) * scale; 222 vector3 y_vec (x,y,z); 223 224 ifs_cont.getline(buffer,BUFF_SIZE); // Z_Vec vector 225 tokenize(vs, buffer); 226 x = atof(vs.at(0).c_str()) * scale; 227 y = atof(vs.at(1).c_str()) * scale; 228 z = atof(vs.at(2).c_str()) * scale; 229 vector3 z_vec (x,y,z); 230 231 // Build unit cell 232 OBUnitCell *cell = new OBUnitCell; 233 cell->SetData(x_vec, y_vec, z_vec); 234 cell->SetSpaceGroup(1); 235 pmol->SetData(cell); 236 237 // Next comes either a list of numbers that represent the stoichiometry of 238 // the cell. The numbers are the atom counts for each type, in the order 239 // listed in the POTCAR file. Since VASP 5.2, a line with a list of atomic 240 // element symbols may precede the atom counts. This will be used if the 241 // POTCAR file is not present. If both are present, the data in the POSCAR 242 // or CONTCAR file will be used. 243 ifs_cont.getline(buffer,BUFF_SIZE); 244 tokenize(vs, buffer); 245 bool symbolsInGeometryFile = false; 246 if (vs.size() != 0) { 247 if (vs.at(0).size() != 0) { 248 if (isalpha(static_cast<int>(vs.at(0).at(0))) != 0) { 249 symbolsInGeometryFile = true; 250 } 251 } 252 } 253 254 // If no element data is present anywhere 255 if (needSymbolsInGeometryFile && !symbolsInGeometryFile && 256 // and there are atoms specified 257 vs.size() != 0) { 258 // Abort read 259 pmol->EndModify(); 260 return false; 261 } 262 263 if (symbolsInGeometryFile) { 264 atomTypes.clear(); 265 for (size_t i = 0; i < vs.size(); ++i) { 266 atomTypes.push_back(OpenBabel::OBElements::GetAtomicNum(vs.at(i).c_str())); 267 } 268 // Fetch next line to get stoichiometry 269 ifs_cont.getline(buffer,BUFF_SIZE); 270 tokenize(vs, buffer); 271 } 272 else if (ifs_pot) { 273 // Get atom types from POTCAR 274 while (ifs_pot.getline(buffer,BUFF_SIZE)) { 275 if (strstr(buffer,"VRHFIN")) { 276 str = buffer; 277 size_t start = str.find("=") + 1; 278 size_t end = str.find(":"); 279 str = str.substr(start, end - start); 280 // Clean up whitespace: 281 for (unsigned int i = 0; i < str.size(); i++) 282 if (str.at(i) == ' ') { 283 str.erase(i,1); 284 --i; 285 } 286 atomTypes.push_back(OpenBabel::OBElements::GetAtomicNum(str.c_str())); 287 } 288 } 289 ifs_pot.close(); 290 } 291 292 // Extract and sum the atom counts. The sum is used to parse the atomic 293 // coordinates 294 totalAtoms = 0; 295 for (unsigned int i = 0; i < vs.size(); i++) { 296 int currentCount = atoi(vs.at(i).c_str()); 297 numAtoms.push_back(currentCount); 298 totalAtoms += currentCount; 299 } 300 301 // Do the number of atom types match the number of atom counts? 302 if (numAtoms.size() != atomTypes.size()) { 303 // If not, abort read 304 pmol->EndModify(); 305 return false; 306 } 307 308 // Cartesian or fractional? 309 ifs_cont.getline(buffer,BUFF_SIZE); 310 selective = false; 311 // Set the variable selective accordingly 312 if (buffer[0] == 'S' || buffer[0] == 's') { 313 selective = true; 314 ifs_cont.getline(buffer,BUFF_SIZE); 315 } 316 // [C|c|K|k] indicates cartesian coordinates, anything else (including 317 // an empty string, buffer[0] == 0) indicates fractional coordinates 318 if ( buffer[0] == 'C' || buffer[0] == 'c' || 319 buffer[0] == 'K' || buffer[0] == 'k' ) { 320 cartesian = true; 321 } 322 else { 323 cartesian = false; 324 } 325 326 atomCount = 0; 327 for (unsigned int i = 0; i < totalAtoms; i++) { 328 // Things get a little nasty here. VASP just prints all the coordinates with no 329 // identifying information one after another here. So in the above sections we've 330 // parsed out which atom types and how many of each are present in atomTypes and 331 // numAtoms, respectively. The counters atomIndex and atomCount have the following 332 // roles: atomIndex keeps track of where we are in atomTypes so that we know the 333 // atomic species we're inserting. atomCount tracks how many of the current 334 // atomTypes.at(atomIndex) species have been inserted, so that when we reach 335 // (atomCount >= numAtoms.at(atomIndex) ) we should stop. Phew. 336 ifs_cont.getline(buffer,BUFF_SIZE); // atom location 337 338 // Let's first figure out the atomic number we are dealing with: 339 while (atomCount >= numAtoms.at(atomIndex)) { 340 atomIndex++; 341 atomCount = 0; 342 } 343 344 // If we made it past that check, we have atomic number = atomTypes.at(atomIndex) 345 // Parse the buffer now. 346 tokenize(vs, buffer); 347 atom = pmol->NewAtom(); 348 atom->SetAtomicNum(atomTypes.at(atomIndex)); 349 x = atof((char*)vs[0].c_str()); 350 y = atof((char*)vs[1].c_str()); 351 z = atof((char*)vs[2].c_str()); 352 vector3 coords (x,y,z); 353 if (!cartesian) 354 coords = cell->FractionalToCartesian( coords ); 355 // If we have Cartesian coordinates, we need to apply the scaling factor 356 else coords *= scale; 357 atom->SetVector(coords); 358 //if the selective dynamics info is present then read it into OBPairData 359 //this needs to be kept somehow to be able to write out the same as input 360 //it's string so it wastes memory :( 361 if (selective && vs.size() >= 6) { 362 key = "move"; 363 value = " "; value += vs[3].c_str(); 364 value += " "; value += vs[4].c_str(); 365 value += " "; value += vs[5].c_str(); 366 cp = new OBPairData; 367 cp->SetAttribute(key); 368 cp->SetValue(value); 369 cp->SetOrigin(fileformatInput); 370 atom->SetData(cp); 371 } 372 373 atomCount++; 374 }; 375 376 // There is some trailing garbage, but AFAIK it's not useful for anything. 377 ifs_cont.close(); 378 379 // Read density of states info from DOSCAR, if available 380 if (ifs_dos) { 381 // Create DOS object 382 OBDOSData *dos = new OBDOSData(); 383 384 // skip header 385 ifs_dos.getline(buffer,BUFF_SIZE); // Junk 386 ifs_dos.getline(buffer,BUFF_SIZE); // Junk 387 ifs_dos.getline(buffer,BUFF_SIZE); // Junk 388 ifs_dos.getline(buffer,BUFF_SIZE); // Junk 389 ifs_dos.getline(buffer,BUFF_SIZE); // Junk 390 391 // Get fermi level 392 double fermi; 393 if (ifs_dos.getline(buffer,BUFF_SIZE)) { // startE endE res fermi ??? 394 tokenize(vs, buffer); 395 fermi = atof(vs[3].c_str()); 396 } 397 398 // Start pulling out energies and densities 399 std::vector<double> energies; 400 std::vector<double> densities; 401 std::vector<double> integration; 402 while (ifs_dos.getline(buffer,BUFF_SIZE)) { 403 tokenize(vs, buffer); 404 energies.push_back(atof(vs[0].c_str())); 405 densities.push_back(atof(vs[1].c_str())); 406 integration.push_back(atof(vs[2].c_str())); 407 } 408 409 if (energies.size() != 0) { 410 dos->SetData(fermi, energies, densities, integration); 411 pmol->SetData(dos); 412 } 413 } 414 415 ifs_dos.close(); 416 417 // Vibration intensities 418 vector3 prevDm; 419 vector<vector3> prevXyz; 420 vector3 currDm; 421 vector<vector3> currXyz; 422 423 // Read in optional information from outcar 424 if (ifs_out) { 425 while (ifs_out.getline(buffer,BUFF_SIZE)) { 426 // Enthalphy 427 if (strstr(buffer, "enthalpy is")) { 428 hasEnthalpy = true; 429 tokenize(vs, buffer); 430 enthalpy_eV = atof(vs[4].c_str()); 431 pv_eV = atof(vs[8].c_str()); 432 } 433 434 // Free energy 435 if (strstr(buffer, "free energy")) { 436 tokenize(vs, buffer); 437 pmol->SetEnergy(atof(vs[4].c_str()) * EV_TO_KCAL_PER_MOL); 438 } 439 440 // Frequencies 441 if (strstr(buffer, "Eigenvectors") && Frequencies.size() == 0) { 442 hasVibrations = true; 443 double x, y, z; 444 ifs_out.getline(buffer,BUFF_SIZE); // dash line 445 ifs_out.getline(buffer,BUFF_SIZE); // blank line 446 ifs_out.getline(buffer,BUFF_SIZE); // blank line 447 ifs_out.getline(buffer,BUFF_SIZE); // first frequency line 448 while (!strstr(buffer, "Eigenvectors")) { 449 vector<vector3> vib; 450 tokenize(vs, buffer); 451 int freqnum = atoi(vs[0].c_str()); 452 if (vs[1].size() == 1 and vs[1].compare("f") == 0) { 453 // Real frequency 454 Frequencies.push_back(atof(vs[7].c_str())); 455 } else if (strstr(vs[1].c_str(), "f/i=")) { 456 // Imaginary frequency 457 Frequencies.push_back(-atof(vs[6].c_str())); 458 } else { 459 // No more frequencies 460 break; 461 } 462 ifs_out.getline(buffer,BUFF_SIZE); // header line 463 ifs_out.getline(buffer,BUFF_SIZE); // first displacement line 464 tokenize(vs, buffer); 465 // normal modes 466 while (vs.size() == 6) { 467 x = atof(vs[3].c_str()); 468 y = atof(vs[4].c_str()); 469 z = atof(vs[5].c_str()); 470 vib.push_back(vector3(x, y, z)); 471 ifs_out.getline(buffer,BUFF_SIZE); // next displacement line 472 tokenize(vs, buffer); 473 } 474 Lx.push_back(vib); 475 ifs_out.getline(buffer,BUFF_SIZE); // next frequency line 476 } 477 } 478 479 if (strstr(buffer, "dipolmoment")) { 480 tokenize(vs, buffer); 481 x = atof(vs[1].c_str()); 482 y = atof(vs[2].c_str()); 483 z = atof(vs[3].c_str()); 484 currDm.Set(x, y, z); 485 } 486 if (strstr(buffer, "TOTAL-FORCE")) { 487 currXyz.clear(); 488 ifs_out.getline(buffer, BUFF_SIZE); // header line 489 ifs_out.getline(buffer, BUFF_SIZE); 490 tokenize(vs, buffer); 491 while (vs.size() == 6) { 492 x = atof(vs[0].c_str()); 493 y = atof(vs[1].c_str()); 494 z = atof(vs[2].c_str()); 495 currXyz.push_back(vector3(x, y, z)); 496 ifs_out.getline(buffer, BUFF_SIZE); // next line 497 tokenize(vs, buffer); 498 } 499 } 500 if (strstr(buffer, "BORN EFFECTIVE CHARGES")) { 501 // IBRION = 7; IBRION = 8 502 dipGrad.clear(); 503 ifs_out.getline(buffer, BUFF_SIZE); // header line 504 ifs_out.getline(buffer, BUFF_SIZE); // `ion #` 505 tokenize(vs, buffer); 506 while (vs.size() == 2) { 507 matrix3x3 dmudq; 508 for (int row = 0; row < 3; ++row) { 509 ifs_out.getline(buffer, BUFF_SIZE); 510 tokenize(vs, buffer); 511 x = atof(vs[1].c_str()); 512 y = atof(vs[2].c_str()); 513 z = atof(vs[3].c_str()); 514 dmudq.SetRow(row, vector3(x, y, z)); 515 } 516 dipGrad.push_back(dmudq); 517 ifs_out.getline(buffer, BUFF_SIZE); // next line 518 tokenize(vs, buffer); 519 } 520 } else if (strstr(buffer, "free energy")) { 521 // IBRION = 5 522 // reached the end of an iteration, use the values 523 if (dipGrad.empty()) { 524 // first iteration: nondisplaced ions 525 dipGrad.resize(pmol->NumAtoms()); 526 } else if (prevXyz.empty()) { 527 // even iteration: store values 528 prevXyz = currXyz; 529 prevDm = currDm; 530 } else { 531 // odd iteration: compute dipGrad = dmu / dxyz for moved ion 532 for (size_t natom = 0; natom < pmol->NumAtoms(); ++natom) { 533 const vector3 dxyz = currXyz[natom] - prevXyz[natom]; 534 vector3::const_iterator iter = std::find_if(dxyz.begin(), dxyz.end(), 535 std::bind2nd(std::not_equal_to<double>(), 0.0)); 536 if (iter != dxyz.end()) dipGrad[natom].SetRow(iter - dxyz.begin(), 537 (currDm - prevDm) / *iter); 538 } 539 prevXyz.clear(); 540 } 541 } 542 } 543 } 544 ifs_out.close(); 545 546 // Set enthalpy 547 if (hasEnthalpy) { 548 OBPairData *enthalpyPD = new OBPairData(); 549 OBPairData *enthalpyPD_pv = new OBPairData(); 550 OBPairData *enthalpyPD_eV = new OBPairData(); 551 OBPairData *enthalpyPD_pv_eV = new OBPairData(); 552 enthalpyPD->SetAttribute("Enthalpy (kcal/mol)"); 553 enthalpyPD_pv->SetAttribute("Enthalpy PV term (kcal/mol)"); 554 enthalpyPD_eV->SetAttribute("Enthalpy (eV)"); 555 enthalpyPD_pv_eV->SetAttribute("Enthalpy PV term (eV)"); 556 double en_kcal_per_mole = enthalpy_eV * EV_TO_KCAL_PER_MOL; 557 double pv_kcal_per_mole = pv_eV * EV_TO_KCAL_PER_MOL; 558 snprintf(tag, BUFF_SIZE, "%f", en_kcal_per_mole); 559 enthalpyPD->SetValue(tag); 560 snprintf(tag, BUFF_SIZE, "%f", pv_kcal_per_mole); 561 enthalpyPD_pv->SetValue(tag); 562 snprintf(tag, BUFF_SIZE, "%f", enthalpy_eV); 563 enthalpyPD_eV->SetValue(tag); 564 snprintf(tag, BUFF_SIZE, "%f", pv_eV); 565 enthalpyPD_pv_eV->SetValue(tag); 566 pmol->SetData(enthalpyPD); 567 pmol->SetData(enthalpyPD_pv); 568 pmol->SetData(enthalpyPD_eV); 569 pmol->SetData(enthalpyPD_pv_eV); 570 } 571 572 // Set vibrations 573 if (hasVibrations) { 574 // compute dDip/dQ 575 vector<double> Intensities; 576 for (vector<vector<vector3> >::const_iterator 577 lxIter = Lx.begin(); lxIter != Lx.end(); ++lxIter) { 578 vector3 intensity; 579 for (size_t natom = 0; natom < dipGrad.size(); ++natom) { 580 intensity += dipGrad[natom].transpose() * lxIter->at(natom) 581 / sqrt(pmol->GetAtomById(natom)->GetAtomicMass()); 582 } 583 Intensities.push_back(dot(intensity, intensity)); 584 } 585 const double max = *max_element(Intensities.begin(), Intensities.end()); 586 if (max != 0.0) { 587 // Normalize 588 std::transform(Intensities.begin(), Intensities.end(), Intensities.begin(), 589 std::bind2nd(std::divides<double>(), max / 100.0)); 590 } else { 591 Intensities.clear(); 592 } 593 OBVibrationData* vd = new OBVibrationData; 594 vd->SetData(Lx, Frequencies, Intensities); 595 pmol->SetData(vd); 596 } 597 598 pmol->EndModify(); 599 600 const char *noBonding = pConv->IsOption("b", OBConversion::INOPTIONS); 601 const char *singleOnly = pConv->IsOption("s", OBConversion::INOPTIONS); 602 603 if (noBonding == nullptr) { 604 pmol->ConnectTheDots(); 605 if (singleOnly == nullptr) { 606 pmol->PerceiveBondOrders(); 607 } 608 } 609 610 return true; 611 } 612 WriteMolecule(OBBase * pOb,OBConversion * pConv)613 bool VASPFormat::WriteMolecule(OBBase* pOb, OBConversion* pConv) 614 { 615 //No surprises in this routine, cartesian coordinates are written out 616 //and if at least a single atom has information about constraints, 617 //then selective dynamics is used and the info is written out. 618 //The atoms are ordered according to their atomic number so that the 619 //output looks nice, this can be reversed by using command line flag "-xw". 620 // 621 OBMol* pmol = dynamic_cast<OBMol*>(pOb); 622 if (pmol == nullptr) { 623 return false; 624 } 625 626 ostream& ofs = *pConv->GetOutStream(); 627 OBMol &mol = *pmol; 628 629 char buffer[BUFF_SIZE]; 630 OBUnitCell *uc = nullptr; 631 vector<vector3> cell; 632 633 const char * sortAtomsNum = pConv->IsOption("w", OBConversion::OUTOPTIONS); 634 const char * sortAtomsCustom = pConv->IsOption("z", OBConversion::OUTOPTIONS); 635 636 // Create a list of ids. These may be sorted by atomic number depending 637 // on the value of keepOrder. 638 std::vector<OBAtom *> atoms_sorted; 639 atoms_sorted.reserve(mol.NumAtoms()); 640 641 FOR_ATOMS_OF_MOL(atom, mol) { 642 atoms_sorted.push_back(&(*atom)); 643 } 644 645 std::vector<int> custom_sort_nums; 646 647 if (sortAtomsCustom != nullptr) 648 { 649 vector<string> vs; 650 tokenize(vs, sortAtomsCustom); 651 for(size_t i = 0; i < vs.size(); ++i) 652 custom_sort_nums.push_back(OBElements::GetAtomicNum(vs[i].c_str())); 653 } 654 655 compare_sort_items csi(custom_sort_nums, sortAtomsNum != nullptr); 656 std::stable_sort(atoms_sorted.begin(), atoms_sorted.end(), csi); 657 658 // Use the atomicNums vector to determine the composition line. 659 // atomicNumsCondensed and atomCounts contain the same data as atomicNums: 660 // if: 661 // atoms_sorted[i]->GetAtomicNum() = [ 3 3 3 2 2 8 2 6 6 ] 662 // then: 663 // atomicNums = [(3 3) (2 2) (8 1) (2 1) (6 2)] 664 665 std::vector<std::pair<int, int> > atomicNums; 666 667 int prev_anum = -20; //not a periodic table number 668 for(int i = 0; i < atoms_sorted.size(); i++) 669 { 670 const int anum = atoms_sorted[i]->GetAtomicNum(); 671 672 if( prev_anum != anum ) 673 { 674 std::pair<int, int> x(anum, 1); 675 atomicNums.push_back(x); 676 } 677 else 678 { 679 if(atomicNums.size() > 0) 680 atomicNums.rbegin()->second++; 681 } 682 683 prev_anum = anum; 684 } 685 686 // write title 687 ofs << mol.GetTitle() << endl; 688 // write the multiplication factor, set this to one 689 // and write the cell using the 3x3 cell matrix 690 ofs << "1.000 " << endl; 691 692 if (!mol.HasData(OBGenericDataType::UnitCell)) { 693 // the unit cell has not been defined. Leave as all zeros so the user 694 // can fill it in themselves 695 for (int ii = 0; ii < 3; ii++) { 696 snprintf(buffer, BUFF_SIZE, "0.0 0.0 0.0"); 697 ofs << buffer << endl; 698 } 699 } 700 else 701 { 702 // there is a unit cell, write it out 703 uc = static_cast<OBUnitCell*>(mol.GetData(OBGenericDataType::UnitCell)); 704 cell = uc->GetCellVectors(); 705 for (vector<vector3>::const_iterator i = cell.begin(); 706 i != cell.end(); ++i) { 707 snprintf(buffer, BUFF_SIZE, "%20.15f%20.15f%20.15f", 708 i->x(), i->y(), i->z()); 709 ofs << buffer << endl; 710 } 711 } 712 713 // go through the atoms first to write out the element names if using 714 // VASP 5 format 715 const char *vasp4Format = pConv->IsOption("4", OBConversion::OUTOPTIONS); 716 if (!vasp4Format) { 717 for (vector< std::pair<int, int> >::const_iterator 718 it = atomicNums.begin(), 719 it_end = atomicNums.end(); it != it_end; ++it) { 720 snprintf(buffer, BUFF_SIZE, "%-3s ", OBElements::GetSymbol(it->first)); 721 ofs << buffer ; 722 } 723 ofs << endl; 724 } 725 726 // then do the same to write out the number of ions of each element 727 for (vector< std::pair<int, int> >::const_iterator 728 it = atomicNums.begin(), 729 it_end = atomicNums.end(); it != it_end; ++it) { 730 snprintf(buffer, BUFF_SIZE, "%-3u ", it->second); 731 ofs << buffer ; 732 } 733 ofs << endl; 734 735 // assume that there are no constraints on the atoms 736 bool selective = false; 737 // and test if any of the atoms has constraints 738 FOR_ATOMS_OF_MOL(atom, mol) { 739 if (atom->HasData("move")){ 740 selective = true; 741 break; 742 } 743 } 744 if (selective) { 745 ofs << "SelectiveDyn" << endl; 746 } 747 748 // print the atomic coordinates in \AA 749 ofs << "Cartesian" << endl; 750 751 for (std::vector<OBAtom *>::const_iterator it = atoms_sorted.begin(); 752 it != atoms_sorted.end(); ++it) 753 { 754 // Print coordinates 755 snprintf(buffer,BUFF_SIZE, "%26.19f %26.19f %26.19f", 756 (*it)->GetX(), (*it)->GetY(), (*it)->GetZ()); 757 ofs << buffer; 758 759 // if at least one atom has info about constraints 760 if (selective) { 761 // if this guy has, write it out 762 if ((*it)->HasData("move")) { 763 OBPairData *cp = (OBPairData*)(*it)->GetData("move"); 764 // seemingly ridiculous number of digits is written out 765 // but sometimes you just don't want to change them 766 ofs << " " << cp->GetValue().c_str(); 767 } 768 else { 769 // the atom has been created and the info has not been copied 770 ofs << " T T T"; 771 } 772 } 773 ofs << endl; 774 } 775 776 return true; 777 } 778 779 } //namespace OpenBabel 780 781