1
2 /*
3 * Copyright 2009-2020 The VOTCA Development Team
4 * (http://www.votca.org)
5 *
6 * Licensed under the Apache License, Version 2.0 (the "License")
7 *
8 * You may not use this file except in compliance with the License.
9 * You may obtain a copy of the License at
10 *
11 * http://www.apache.org/licenses/LICENSE-2.0
12 *
13 * Unless required by applicable law or agreed to in writing, software
14 * distributed under the License is distributed on an "AS IS" BASIS,
15 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
16 * See the License for the specific language governing permissions and
17 * limitations under the License.
18 *
19 */
20
21 // Standard includes
22 #include <cstdio>
23 #include <iomanip>
24
25 // Third party includes
26 #include <boost/algorithm/string.hpp>
27 #include <boost/filesystem.hpp>
28 #include <boost/format.hpp>
29
30 // VOTCA includes
31 #include <stdexcept>
32 #include <string>
33 #include <votca/tools/elements.h>
34 #include <votca/tools/getline.h>
35
36 // Local VOTCA includes
37 #include "votca/tools/globals.h"
38 #include "votca/tools/tokenizer.h"
39 #include "votca/xtp/basisset.h"
40 #include "votca/xtp/ecpaobasis.h"
41 #include "votca/xtp/molden.h"
42 #include "votca/xtp/orbitals.h"
43
44 // Local private VOTCA includes
45 #include "orca.h"
46
47 namespace votca {
48 namespace xtp {
49 using namespace std;
50
ParseSpecificOptions(const tools::Property & options)51 void Orca::ParseSpecificOptions(const tools::Property& options) {
52
53 // Orca file names
54 std::string fileName = options.get("temporary_file").as<std::string>();
55
56 input_file_name_ = fileName + ".inp";
57 log_file_name_ = fileName + ".log";
58 shell_file_name_ = fileName + ".sh";
59 mo_file_name_ = fileName + ".gbw";
60 }
61
62 /* Custom basis sets are written on a per-element basis to
63 * the system.bas/aux file(s), which are then included in the
64 * Orca input file using GTOName = "system.bas/aux"
65 */
WriteBasisset(const QMMolecule & qmatoms,std::string & bs_name,std::string & el_file_name)66 void Orca::WriteBasisset(const QMMolecule& qmatoms, std::string& bs_name,
67 std::string& el_file_name) {
68
69 std::vector<std::string> UniqueElements = qmatoms.FindUniqueElements();
70
71 tools::Elements elementInfo;
72 BasisSet bs;
73 bs.Load(bs_name);
74 XTP_LOG(Log::error, *pLog_) << "Loaded Basis Set " << bs_name << flush;
75 ofstream el_file;
76
77 el_file.open(el_file_name);
78 el_file << "$DATA" << endl;
79
80 for (const std::string& element_name : UniqueElements) {
81 const Element& element = bs.getElement(element_name);
82 el_file << elementInfo.getEleFull(element_name) << endl;
83 for (const Shell& shell : element) {
84 el_file << xtp::EnumToString(shell.getL()) << " " << shell.getSize()
85 << endl;
86 Index sh_idx = 0;
87 for (const GaussianPrimitive& gaussian : shell) {
88 sh_idx++;
89 el_file << " " << sh_idx << " " << indent(gaussian.decay());
90 el_file << " " << indent(gaussian.contraction());
91 el_file << endl;
92 }
93 }
94 }
95 el_file << "STOP\n";
96 el_file.close();
97
98 return;
99 }
100
101 /* Coordinates are written in standard Element,x,y,z format to the
102 * input file.
103 */
WriteCoordinates(std::ofstream & inp_file,const QMMolecule & qmatoms)104 void Orca::WriteCoordinates(std::ofstream& inp_file,
105 const QMMolecule& qmatoms) {
106
107 for (const QMAtom& atom : qmatoms) {
108 Eigen::Vector3d pos = atom.getPos() * tools::conv::bohr2ang;
109 inp_file << setw(3) << atom.getElement() << setw(12)
110 << setiosflags(ios::fixed) << setprecision(6) << pos.x()
111 << setw(12) << setiosflags(ios::fixed) << setprecision(6)
112 << pos.y() << setw(12) << setiosflags(ios::fixed)
113 << setprecision(6) << pos.z() << endl;
114 }
115 inp_file << "* \n" << endl;
116 return;
117 }
118
119 /* If custom ECPs are used, they need to be specified in the input file
120 * in a section following the basis set includes.
121 */
WriteECP(std::ofstream & inp_file,const QMMolecule & qmatoms)122 void Orca::WriteECP(std::ofstream& inp_file, const QMMolecule& qmatoms) {
123
124 inp_file << endl;
125 std::vector<std::string> UniqueElements = qmatoms.FindUniqueElements();
126
127 ECPBasisSet ecp;
128 ecp.Load(options_.get("ecp").as<std::string>());
129
130 XTP_LOG(Log::error, *pLog_) << "Loaded Pseudopotentials "
131 << options_.get("ecp").as<std::string>() << flush;
132
133 for (const std::string& element_name : UniqueElements) {
134 try {
135 ecp.getElement(element_name);
136 } catch (std::runtime_error& error) {
137 XTP_LOG(Log::error, *pLog_)
138 << "No pseudopotential for " << element_name << " available" << flush;
139 continue;
140 }
141 const ECPElement& element = ecp.getElement(element_name);
142
143 inp_file << "\n"
144 << "NewECP"
145 << " " << element_name << endl;
146 inp_file << "N_core"
147 << " " << element.getNcore() << endl;
148 inp_file << "lmax"
149 << " " << EnumToString(element.getLmax()) << endl;
150 // For Orca the order doesn't matter but let's write it in ascending order
151 // write remaining shells in ascending order s,p,d...
152 for (Index i = 0; i <= Index(element.getLmax()); i++) {
153 for (const ECPShell& shell : element) {
154 if (Index(shell.getL()) == i) {
155 // shell type, number primitives, scale factor
156 inp_file << xtp::EnumToString(shell.getL()) << " " << shell.getSize()
157 << endl;
158 Index sh_idx = 0;
159 for (const ECPGaussianPrimitive& gaussian : shell) {
160 sh_idx++;
161 inp_file << sh_idx << " " << gaussian.decay_ << " "
162 << gaussian.contraction_ << " " << gaussian.power_ << endl;
163 }
164 }
165 }
166 }
167 inp_file << "end\n "
168 << "\n"
169 << endl;
170 }
171 return;
172 }
173
WriteChargeOption()174 void Orca::WriteChargeOption() {
175 this->options_.addTree("orca.pointcharges", "\"background.crg\"");
176 }
177
178 /* For QM/MM the molecules in the MM environment are represented by
179 * their atomic partial charge distributions. ORCA expects them in
180 * q,x,y,z format in a separate file "background.crg"
181 */
WriteBackgroundCharges()182 void Orca::WriteBackgroundCharges() {
183
184 std::ofstream crg_file;
185 std::string crg_file_name_full_ = run_dir_ + "/background.crg";
186 crg_file.open(crg_file_name_full_);
187 Index total_background = 0;
188
189 for (const std::unique_ptr<StaticSite>& site : externalsites_) {
190 if (site->getCharge() != 0.0) {
191 total_background++;
192 }
193 total_background += SplitMultipoles(*site).size();
194 } // counting only
195
196 crg_file << total_background << endl;
197 boost::format fmt("%1$+1.7f %2$+1.7f %3$+1.7f %4$+1.7f");
198 // now write
199 for (const std::unique_ptr<StaticSite>& site : externalsites_) {
200 Eigen::Vector3d pos = site->getPos() * tools::conv::bohr2ang;
201 string sitestring =
202 boost::str(fmt % site->getCharge() % pos.x() % pos.y() % pos.z());
203 if (site->getCharge() != 0.0) {
204 crg_file << sitestring << endl;
205 }
206 std::vector<MinimalMMCharge> split_multipoles = SplitMultipoles(*site);
207 for (const auto& mpoles : split_multipoles) {
208 Eigen::Vector3d pos2 = mpoles.pos_ * tools::conv::bohr2ang;
209 string multipole =
210 boost::str(fmt % mpoles.q_ % pos2.x() % pos2.y() % pos2.z());
211 crg_file << multipole << endl;
212 }
213 }
214
215 return;
216 }
217
218 /**
219 * Prepares the *.inp file from a vector of segments
220 * Appends a guess constructed from monomer orbitals if supplied, Not
221 * implemented yet
222 */
WriteInputFile(const Orbitals & orbitals)223 bool Orca::WriteInputFile(const Orbitals& orbitals) {
224
225 std::vector<std::string> results;
226 std::string temp_suffix = "/id";
227 std::string scratch_dir_backup = scratch_dir_;
228 std::ofstream inp_file;
229 std::string inp_file_name_full = run_dir_ + "/" + input_file_name_;
230 inp_file.open(inp_file_name_full);
231 // header
232 inp_file << "* xyz " << charge_ << " " << spin_ << endl;
233 Index threads = OPENMP::getMaxThreads();
234 const QMMolecule& qmatoms = orbitals.QMAtoms();
235 // put coordinates
236 WriteCoordinates(inp_file, qmatoms);
237 // add parallelization info
238 inp_file << "%pal\n"
239 << "nprocs " << threads << "\nend"
240 << "\n"
241 << endl;
242 // basis set info
243 std::string el_file_name = run_dir_ + "/" + "system.bas";
244 WriteBasisset(qmatoms, basisset_name_, el_file_name);
245 inp_file << "%basis\n";
246 inp_file << "GTOName"
247 << " "
248 << "="
249 << "\"system.bas\";" << endl;
250 if (options_.exists("auxbasisset")) {
251 std::string aux_file_name = run_dir_ + "/" + "system.aux";
252 std::string auxbasisset_name =
253 options_.get("auxbasisset").as<std::string>();
254 WriteBasisset(qmatoms, auxbasisset_name, aux_file_name);
255 inp_file << "GTOAuxName"
256 << " "
257 << "="
258 << "\"system.aux\";" << endl;
259 } // write_auxbasis set
260
261 // ECPs
262 if (options_.exists("ecp")) {
263 WriteECP(inp_file, qmatoms);
264 }
265 inp_file << "end\n "
266 << "\n"
267 << endl; // This end is for the basis set block
268 if (!externalsites_.empty()) {
269 WriteBackgroundCharges();
270 }
271
272 // External Electric field
273 if (options_.exists("externalfield")) {
274 Eigen::Vector3d field = options_.get("externalfield").as<Eigen::Vector3d>();
275 inp_file << "%scf\n ";
276 inp_file << " efield " << field.x() << ", " << field.y() << ", "
277 << field.z() << "\n";
278 inp_file << "end\n";
279 inp_file << std::endl;
280 }
281 std::string input_options;
282 // Write Orca section specified by the user
283 for (const auto& prop : this->options_.get("orca")) {
284 const std::string& prop_name = prop.name();
285 if (prop_name == "pointcharges") {
286 input_options += this->CreateInputSection("orca.pointcharges");
287 } else if (prop_name != "method") {
288 input_options += this->CreateInputSection("orca." + prop_name);
289 }
290 }
291 // Write main DFT method
292 input_options += this->WriteMethod();
293 inp_file << input_options;
294 inp_file << '\n';
295 inp_file.close();
296 // and now generate a shell script to run both jobs, if neccessary
297
298 XTP_LOG(Log::info, *pLog_)
299 << "Setting the scratch dir to " << scratch_dir_ + temp_suffix << flush;
300 scratch_dir_ = scratch_dir_backup + temp_suffix;
301 WriteShellScript();
302 scratch_dir_ = scratch_dir_backup;
303 return true;
304 }
305
WriteShellScript()306 bool Orca::WriteShellScript() {
307 ofstream shell_file;
308 std::string shell_file_name_full = run_dir_ + "/" + shell_file_name_;
309 shell_file.open(shell_file_name_full);
310 shell_file << "#!/bin/bash" << endl;
311 shell_file << "mkdir -p " << scratch_dir_ << endl;
312
313 if (options_.get("initial_guess").as<std::string>() == "orbfile") {
314 if (!(boost::filesystem::exists(run_dir_ + "/molA.gbw") &&
315 boost::filesystem::exists(run_dir_ + "/molB.gbw"))) {
316 throw runtime_error(
317 "Using guess relies on a molA.gbw and a molB.gbw file being in the "
318 "directory.");
319 }
320 shell_file << options_.get("executable").as<std::string>()
321 << "_mergefrag molA.gbw molB.gbw dimer.gbw > merge.log" << endl;
322 }
323 shell_file << options_.get("executable").as<std::string>() << " "
324 << input_file_name_ << " > " << log_file_name_
325 << endl; //" 2> run.error" << endl;
326 std::string base_name = mo_file_name_.substr(0, mo_file_name_.size() - 4);
327 shell_file << options_.get("executable").as<std::string>() << "_2mkl "
328 << base_name << " -molden" << endl;
329 shell_file.close();
330 return true;
331 }
332
333 /**
334 * Runs the Orca job.
335 */
RunDFT()336 bool Orca::RunDFT() {
337
338 XTP_LOG(Log::error, *pLog_) << "Running Orca job\n" << flush;
339
340 if (std::system(nullptr)) {
341
342 std::string command = "cd " + run_dir_ + "; sh " + shell_file_name_;
343 Index check = std::system(command.c_str());
344 if (check == -1) {
345 XTP_LOG(Log::error, *pLog_)
346 << input_file_name_ << " failed to start" << flush;
347 return false;
348 }
349 if (CheckLogFile()) {
350 XTP_LOG(Log::error, *pLog_) << "Finished Orca job" << flush;
351 return true;
352 } else {
353 XTP_LOG(Log::error, *pLog_) << "Orca job failed" << flush;
354 }
355 } else {
356 XTP_LOG(Log::error, *pLog_)
357 << input_file_name_ << " failed to start" << flush;
358 return false;
359 }
360
361 return true;
362 }
363
364 /**
365 * Cleans up after the Orca job
366 */
CleanUp()367 void Orca::CleanUp() {
368
369 if (options_.get("initial_guess").as<std::string>() == "orbfile") {
370 remove((run_dir_ + "/" + "molA.gbw").c_str());
371 remove((run_dir_ + "/" + "molB.gbw").c_str());
372 remove((run_dir_ + "/" + "dimer.gbw").c_str());
373 }
374 // cleaning up the generated files
375 if (!cleanup_.empty()) {
376 tools::Tokenizer tok_cleanup(cleanup_, ",");
377 for (const std::string& substring : tok_cleanup) {
378 if (substring == "inp") {
379 std::string file_name = run_dir_ + "/" + input_file_name_;
380 remove(file_name.c_str());
381 }
382
383 if (substring == "bas") {
384 std::string file_name = run_dir_ + "/system.bas";
385 remove(file_name.c_str());
386 }
387
388 if (substring == "log") {
389 std::string file_name = run_dir_ + "/" + log_file_name_;
390 remove(file_name.c_str());
391 }
392
393 if (substring == "gbw") {
394 std::string file_name = run_dir_ + "/" + mo_file_name_;
395 remove(file_name.c_str());
396 }
397
398 if (substring == "ges") {
399 std::string file_name = run_dir_ + "/system.ges";
400 remove(file_name.c_str());
401 }
402 if (substring == "prop") {
403 std::string file_name = run_dir_ + "/system.prop";
404 remove(file_name.c_str());
405 }
406 }
407 }
408 return;
409 }
410
GetCharges() const411 StaticSegment Orca::GetCharges() const {
412
413 StaticSegment result("charges", 0);
414
415 XTP_LOG(Log::error, *pLog_) << "Parsing " << log_file_name_ << flush;
416 std::string log_file_name_full = run_dir_ + "/" + log_file_name_;
417 std::string line;
418
419 std::ifstream input_file(log_file_name_full);
420 while (input_file) {
421 tools::getline(input_file, line);
422 boost::trim(line);
423 GetCoordinates(result, line, input_file);
424
425 std::string::size_type charge_pos = line.find("CHELPG Charges");
426
427 if (charge_pos != std::string::npos) {
428 XTP_LOG(Log::error, *pLog_) << "Getting charges" << flush;
429 tools::getline(input_file, line);
430 std::vector<std::string> row = GetLineAndSplit(input_file, "\t ");
431 Index nfields = Index(row.size());
432 bool hasAtoms = result.size() > 0;
433 while (nfields == 4) {
434 Index atom_id = boost::lexical_cast<Index>(row.at(0));
435 std::string atom_type = row.at(1);
436 double atom_charge = boost::lexical_cast<double>(row.at(3));
437 row = GetLineAndSplit(input_file, "\t ");
438 nfields = Index(row.size());
439 if (hasAtoms) {
440 StaticSite& temp = result.at(atom_id);
441 if (temp.getElement() != atom_type) {
442 throw std::runtime_error(
443 "Getting charges failed. Mismatch in elemts:" +
444 temp.getElement() + " vs " + atom_type);
445 }
446 temp.setCharge(atom_charge);
447 } else {
448 StaticSite temp =
449 StaticSite(atom_id, atom_type, Eigen::Vector3d::Zero());
450 temp.setCharge(atom_charge);
451 result.push_back(temp);
452 }
453 }
454 }
455 }
456 return result;
457 }
458
GetPolarizability() const459 Eigen::Matrix3d Orca::GetPolarizability() const {
460 std::string line;
461 ifstream input_file((run_dir_ + "/" + log_file_name_));
462 bool has_pol = false;
463
464 Eigen::Matrix3d pol = Eigen::Matrix3d::Zero();
465 while (input_file) {
466 tools::getline(input_file, line);
467 boost::trim(line);
468
469 std::string::size_type pol_pos = line.find("THE POLARIZABILITY TENSOR");
470 if (pol_pos != std::string::npos) {
471 XTP_LOG(Log::error, *pLog_) << "Getting polarizability" << flush;
472 tools::getline(input_file, line);
473 tools::getline(input_file, line);
474 tools::getline(input_file, line);
475
476 if (line.find("The raw cartesian tensor (atomic units)") ==
477 std::string::npos) {
478 throw std::runtime_error(
479 "Could not find cartesian polarization tensor");
480 }
481
482 for (Index i = 0; i < 3; i++) {
483 tools::getline(input_file, line);
484 std::vector<double> values =
485 tools::Tokenizer(line, " ").ToVector<double>();
486 if (values.size() != 3) {
487 throw std::runtime_error("polarization line " + line +
488 " cannot be parsed");
489 }
490 Eigen::Vector3d row;
491 row << values[0], values[1], values[2];
492 pol.row(i) = row;
493 }
494
495 has_pol = true;
496 }
497 }
498 if (!has_pol) {
499 throw std::runtime_error("Could not find polarization in logfile");
500 }
501 return pol;
502 }
503
ParseLogFile(Orbitals & orbitals)504 bool Orca::ParseLogFile(Orbitals& orbitals) {
505 bool found_success = false;
506 orbitals.setQMpackage(getPackageName());
507 orbitals.setDFTbasisName(basisset_name_);
508 if (options_.exists("ecp")) {
509 orbitals.setECPName(options_.get("ecp").as<std::string>());
510 }
511
512 orbitals.setXCGrid("xfine"); // TODO find a better approximation for the
513 // orca grid.
514 orbitals.setXCFunctionalName(options_.get("functional").as<std::string>());
515
516 XTP_LOG(Log::error, *pLog_) << "Parsing " << log_file_name_ << flush;
517 std::string log_file_name_full = run_dir_ + "/" + log_file_name_;
518 // check if LOG file is complete
519 if (!CheckLogFile()) {
520 return false;
521 }
522 std::map<Index, double> energies;
523 std::map<Index, double> occupancy;
524
525 std::string line;
526 Index levels = 0;
527 Index number_of_electrons = 0;
528 std::vector<std::string> results;
529
530 std::ifstream input_file(log_file_name_full);
531
532 if (input_file.fail()) {
533 XTP_LOG(Log::error, *pLog_)
534 << "File " << log_file_name_full << " not found " << flush;
535 return false;
536 } else {
537 XTP_LOG(Log::error, *pLog_)
538 << "Reading Coordinates and occupationnumbers and energies from "
539 << log_file_name_full << flush;
540 }
541 // Coordinates of the final configuration depending on whether it is an
542 // optimization or not
543
544 QMMolecule& mol = orbitals.QMAtoms();
545 while (input_file) {
546 tools::getline(input_file, line);
547 boost::trim(line);
548
549 GetCoordinates(mol, line, input_file);
550
551 std::string::size_type energy_pos = line.find("FINAL SINGLE");
552 if (energy_pos != std::string::npos) {
553 results = tools::Tokenizer(line, " ").ToVector();
554 std::string energy = results[4];
555 boost::trim(energy);
556 orbitals.setQMEnergy(boost::lexical_cast<double>(energy));
557 XTP_LOG(Log::error, *pLog_) << (boost::format("QM energy[Hrt]: %4.6f ") %
558 orbitals.getDFTTotalEnergy())
559 .str()
560 << flush;
561 }
562
563 std::string::size_type HFX_pos = line.find("Fraction HF Exchange ScalHFX");
564
565 if (HFX_pos != std::string::npos) {
566 results = tools::Tokenizer(line, " ").ToVector();
567 double ScaHFX = boost::lexical_cast<double>(results.back());
568
569 orbitals.setScaHFX(ScaHFX);
570 XTP_LOG(Log::error, *pLog_)
571 << "DFT with " << ScaHFX << " of HF exchange!" << flush;
572 }
573
574 std::string::size_type dim_pos = line.find("Basis Dimension");
575 if (dim_pos != std::string::npos) {
576 results = tools::Tokenizer(line, " ").ToVector();
577 std::string dim =
578 results[4]; // The 4th element of results vector is the Basis Dim
579 boost::trim(dim);
580 levels = boost::lexical_cast<Index>(dim);
581 XTP_LOG(Log::info, *pLog_) << "Basis Dimension: " << levels << flush;
582 XTP_LOG(Log::info, *pLog_) << "Energy levels: " << levels << flush;
583 }
584
585 std::string::size_type OE_pos = line.find("ORBITAL ENERGIES");
586 if (OE_pos != std::string::npos) {
587
588 number_of_electrons = 0;
589 tools::getline(input_file, line);
590 tools::getline(input_file, line);
591 tools::getline(input_file, line);
592 if (line.find("E(Eh)") == std::string::npos) {
593 XTP_LOG(Log::error, *pLog_)
594 << "Warning: Orbital Energies not found in log file" << flush;
595 }
596 for (Index i = 0; i < levels; i++) {
597 results = GetLineAndSplit(input_file, " ");
598 std::string no = results[0];
599 boost::trim(no);
600 Index levelnumber = boost::lexical_cast<Index>(no);
601 if (levelnumber != i) {
602 XTP_LOG(Log::error, *pLog_) << "Have a look at the orbital energies "
603 "something weird is going on"
604 << flush;
605 }
606 std::string oc = results[1];
607 boost::trim(oc);
608 double occ = boost::lexical_cast<double>(oc);
609 // We only count alpha electrons, each orbital must be empty or doubly
610 // occupied
611 if (occ == 2 || occ == 1) {
612 number_of_electrons++;
613 occupancy[i] = occ;
614 if (occ == 1) {
615 XTP_LOG(Log::error, *pLog_)
616 << "Watch out! No distinction between alpha and beta "
617 "electrons. Check if occ = 1 is suitable for your "
618 "calculation "
619 << flush;
620 }
621 } else if (occ == 0) {
622 occupancy[i] = occ;
623 } else {
624 throw runtime_error(
625 "Only empty or doubly occupied orbitals are allowed not "
626 "running the right kind of DFT calculation");
627 }
628 std::string e = results[2];
629 boost::trim(e);
630 energies[i] = boost::lexical_cast<double>(e);
631 }
632 }
633
634 std::string::size_type success =
635 line.find("* SUCCESS *");
636 if (success != std::string::npos) {
637 found_success = true;
638 }
639 }
640
641 XTP_LOG(Log::info, *pLog_)
642 << "Alpha electrons: " << number_of_electrons << flush;
643 Index occupied_levels = number_of_electrons;
644 Index unoccupied_levels = levels - occupied_levels;
645 XTP_LOG(Log::info, *pLog_) << "Occupied levels: " << occupied_levels << flush;
646 XTP_LOG(Log::info, *pLog_)
647 << "Unoccupied levels: " << unoccupied_levels << flush;
648
649 /************************************************************/
650
651 // copying information to the orbitals object
652
653 orbitals.setBasisSetSize(levels);
654 orbitals.setNumberOfAlphaElectrons(number_of_electrons);
655 orbitals.setNumberOfOccupiedLevels(occupied_levels);
656
657 // copying energies to a vector
658 orbitals.MOs().eigenvalues().resize(levels);
659 // level_ = 1;
660 for (Index i = 0; i < levels; i++) {
661 orbitals.MOs().eigenvalues()[i] = energies[i];
662 }
663
664 XTP_LOG(Log::error, *pLog_) << "Done reading Log file" << flush;
665
666 return found_success;
667 }
668 template <class T>
GetCoordinates(T & mol,string & line,ifstream & input_file) const669 void Orca::GetCoordinates(T& mol, string& line, ifstream& input_file) const {
670 std::string::size_type coordinates_pos =
671 line.find("CARTESIAN COORDINATES (ANGSTROEM)");
672
673 using Atom = typename T::Atom_Type;
674
675 if (coordinates_pos != std::string::npos) {
676 XTP_LOG(Log::error, *pLog_) << "Getting the coordinates" << flush;
677 bool has_QMAtoms = mol.size() > 0;
678 // three garbage lines
679 tools::getline(input_file, line);
680 // now starts the data in format
681 // id_ type Qnuc x y z
682 vector<string> row = GetLineAndSplit(input_file, "\t ");
683 Index nfields = Index(row.size());
684 Index atom_id = 0;
685 while (nfields == 4) {
686 string atom_type = row.at(0);
687 double x = boost::lexical_cast<double>(row.at(1));
688 double y = boost::lexical_cast<double>(row.at(2));
689 double z = boost::lexical_cast<double>(row.at(3));
690 row = GetLineAndSplit(input_file, "\t ");
691 nfields = Index(row.size());
692 Eigen::Vector3d pos(x, y, z);
693 pos *= tools::conv::ang2bohr;
694 if (has_QMAtoms == false) {
695 mol.push_back(Atom(atom_id, atom_type, pos));
696 } else {
697 Atom& pAtom = mol.at(atom_id);
698 pAtom.setPos(pos);
699 }
700 atom_id++;
701 }
702 }
703 }
704
CheckLogFile()705 bool Orca::CheckLogFile() {
706 // check if the log file exists
707 ifstream input_file(run_dir_ + "/" + log_file_name_);
708
709 if (input_file.fail()) {
710 XTP_LOG(Log::error, *pLog_) << "Orca LOG is not found" << flush;
711 return false;
712 };
713
714 std::string line;
715 while (input_file) {
716 tools::getline(input_file, line);
717 boost::trim(line);
718 std::string::size_type error = line.find("FATAL ERROR ENCOUNTERED");
719 if (error != std::string::npos) {
720 XTP_LOG(Log::error, *pLog_) << "ORCA encountered a fatal error, maybe a "
721 "look in the log file may help."
722 << flush;
723 return false;
724 }
725 error = line.find(
726 "mpirun detected that one or more processes exited with non-zero "
727 "status");
728 if (error != std::string::npos) {
729 XTP_LOG(Log::error, *pLog_)
730 << "ORCA had an mpi problem, maybe your openmpi version is not good."
731 << flush;
732 return false;
733 }
734 }
735 return true;
736 }
737
738 // Parses the molden file from orca and stores data in the Orbitals object
ParseMOsFile(Orbitals & orbitals)739 bool Orca::ParseMOsFile(Orbitals& orbitals) {
740 if (!CheckLogFile()) {
741 return false;
742 }
743 std::vector<double> coefficients;
744 Index basis_size = orbitals.getBasisSetSize();
745 if (basis_size == 0) {
746 throw runtime_error(
747 "Basis size not set, calculator does not parse log file first");
748 }
749
750 XTP_LOG(Log::error, *pLog_) << "Reading Molden file" << flush;
751
752 Molden molden(*pLog_);
753
754 if (orbitals.getDFTbasisName() == "") {
755 throw runtime_error(
756 "Basisset names should be set before reading the molden file.");
757 }
758
759 molden.setBasissetInfo(orbitals.getDFTbasisName());
760 std::string file_name = run_dir_ + "/" +
761 mo_file_name_.substr(0, mo_file_name_.size() - 4) +
762 ".molden.input";
763 XTP_LOG(Log::error, *pLog_) << "Molden file: " << file_name << flush;
764 std::ifstream molden_file(file_name);
765 if (!molden_file.good()) {
766 throw std::runtime_error(
767 "Could not find the molden input file for the MO coefficients.\nIf you "
768 "have run the orca calculation manually or use data from an old\n"
769 "calculation, make sure that besides the .gbw file a .molden.input is\n"
770 "present. If not, convert the .gbw file to a .molden.input file with\n"
771 "the orca_2mkl tool from orca.\nAn example, if you have a benzene.gbw "
772 "file run:\n orca_2mkl benzene -molden\n");
773 }
774
775 molden.parseMoldenFile(file_name, orbitals);
776
777 XTP_LOG(Log::error, *pLog_) << "Done parsing" << flush;
778 return true;
779 }
780
indent(const double & number)781 std::string Orca::indent(const double& number) {
782 std::stringstream ssnumber;
783 if (number >= 0) {
784 ssnumber << " ";
785 } else {
786 ssnumber << " ";
787 }
788 ssnumber << setiosflags(ios::fixed) << setprecision(15) << std::scientific
789 << number;
790 std::string snumber = ssnumber.str();
791 return snumber;
792 }
793
CreateInputSection(const std::string & key) const794 std::string Orca::CreateInputSection(const std::string& key) const {
795 std::stringstream stream;
796 std::string section = key.substr(key.find(".") + 1);
797 stream << "%" << section;
798 if (KeywordIsSingleLine(key)) {
799 stream << " " << options_.get(key).as<std::string>() << "\n";
800 } else {
801 stream << "\n"
802 << options_.get(key).as<std::string>() << "\n"
803 << "end\n";
804 }
805
806 return stream.str();
807 }
808
KeywordIsSingleLine(const std::string & key) const809 bool Orca::KeywordIsSingleLine(const std::string& key) const {
810 tools::Tokenizer values(this->options_.get(key).as<std::string>(), " ");
811 std::vector<std::string> words = values.ToVector();
812 return ((words.size() <= 1) ? true : false);
813 }
814
WriteMethod() const815 std::string Orca::WriteMethod() const {
816 std::stringstream stream;
817 std::string opt = (options_.get("optimize").as<bool>()) ? "Opt" : "";
818 const tools::Property& orca = options_.get("orca");
819 std::string user_method =
820 (orca.exists("method")) ? orca.get("method").as<std::string>() : "";
821 std::string convergence = "";
822 if (!orca.exists("scf")) {
823 convergence = this->convergence_map_.at(
824 options_.get("convergence_tightness").as<std::string>());
825 }
826 stream << "! DFT " << this->GetOrcaFunctionalName() << " " << convergence
827 << " " << opt
828 << " "
829 // additional properties provided by the user
830 << user_method << "\n";
831 return stream.str();
832 }
833
GetOrcaFunctionalName() const834 std::string Orca::GetOrcaFunctionalName() const {
835
836 std::map<std::string, std::string> votca_to_orca;
837
838 votca_to_orca["XC_HYB_GGA_XC_B1LYP"] = "B1LYP";
839 votca_to_orca["XC_HYB_GGA_XC_B3LYP"] = "B3LYP";
840 votca_to_orca["XC_HYB_GGA_XC_PBEH"] = "PBE0";
841 votca_to_orca["XC_GGA_C_PBE"] = "PBE";
842 votca_to_orca["XC_GGA_X_PBE"] = "PBE";
843
844 std::string votca_functional =
845 options_.get("functional").as<std::vector<std::string>>()[0];
846
847 std::string orca_functional;
848 if (votca_to_orca.count(votca_functional)) {
849 orca_functional = votca_to_orca.at(votca_functional);
850 } else if (options_.exists("orca." + votca_functional)) {
851 orca_functional =
852 options_.get("orca." + votca_functional).as<std::string>();
853 } else {
854 throw std::runtime_error(
855 "Cannot translate " + votca_functional +
856 " to orca functional names. Please add a <" + votca_functional +
857 "> to your orca input options with the functional orca should use.");
858 }
859 return orca_functional;
860 }
861
862 } // namespace xtp
863 } // namespace votca
864