1 /*
2 * Copyright 2009-2020 The VOTCA Development Team
3 * (http://www.votca.org)
4 *
5 * Licensed under the Apache License, Version 2.0 (the "License")
6 *
7 * You may not use this file except in compliance with the License.
8 * You may obtain a copy of the License at
9 *
10 * http://www.apache.org/licenses/LICENSE-2.0
11 *
12 * Unless required by applicable law or agreed to in writing, software
13 * distributed under the License is distributed on an "AS IS" BASIS,
14 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15 * See the License for the specific language governing permissions and
16 * limitations under the License.
17 *
18 */
19
20 // Third party includes
21 #include <boost/algorithm/string.hpp>
22
23 // Local VOTCA includes
24 #include "votca/tools/getline.h"
25 #include "votca/tools/globals.h"
26 #include "votca/xtp/ecpaobasis.h"
27 #include "votca/xtp/orbitals.h"
28 #include "votca/xtp/qmpackage.h"
29 #include "votca/xtp/qmpackagefactory.h"
30
31 namespace votca {
32 namespace xtp {
33 using std::flush;
34
Initialize(const tools::Property & options)35 void QMPackage::Initialize(const tools::Property& options) {
36 charge_ = options.get("charge").as<Index>();
37 spin_ = options.get("spin").as<Index>();
38 basisset_name_ = options.get("basisset").as<std::string>();
39
40 cleanup_ = options.get("cleanup").as<std::string>();
41 scratch_dir_ = options.get("scratch").as<std::string>();
42
43 options_ = options;
44
45 ParseSpecificOptions(options);
46 }
47
Run()48 bool QMPackage::Run() {
49 std::chrono::time_point<std::chrono::system_clock> start =
50 std::chrono::system_clock::now();
51
52 bool error_value = RunDFT();
53
54 std::chrono::duration<double> elapsed_time =
55 std::chrono::system_clock::now() - start;
56 XTP_LOG(Log::error, *pLog_) << TimeStamp() << " DFT calculation took "
57 << elapsed_time.count() << " seconds." << flush;
58 return error_value;
59 }
60
ReorderOutput(Orbitals & orbitals) const61 void QMPackage::ReorderOutput(Orbitals& orbitals) const {
62 if (!orbitals.hasQMAtoms()) {
63 throw std::runtime_error("Orbitals object has no QMAtoms");
64 }
65
66 AOBasis dftbasis = orbitals.SetupDftBasis();
67 // necessary to update nuclear charges on qmatoms
68 if (orbitals.hasECPName()) {
69 ECPBasisSet ecps;
70 ecps.Load(orbitals.getECPName());
71 ECPAOBasis ecpbasis;
72 ecpbasis.Fill(ecps, orbitals.QMAtoms());
73 }
74
75 if (orbitals.hasMOs()) {
76 OrbReorder reorder(ShellReorder(), ShellMulitplier());
77 reorder.reorderOrbitals(orbitals.MOs().eigenvectors(), dftbasis);
78 XTP_LOG(Log::info, *pLog_) << "Reordered MOs" << flush;
79 }
80
81 return;
82 }
83
ReorderMOsBack(const Orbitals & orbitals) const84 Eigen::MatrixXd QMPackage::ReorderMOsBack(const Orbitals& orbitals) const {
85 if (!orbitals.hasQMAtoms()) {
86 throw std::runtime_error("Orbitals object has no QMAtoms");
87 }
88 AOBasis dftbasis = orbitals.SetupDftBasis();
89 Eigen::MatrixXd result = orbitals.MOs().eigenvectors();
90 bool reverseOrder = true;
91 OrbReorder reorder(ShellReorder(), ShellMulitplier(), reverseOrder);
92 reorder.reorderOrbitals(result, dftbasis);
93 return result;
94 }
95
SplitMultipoles(const StaticSite & aps) const96 std::vector<QMPackage::MinimalMMCharge> QMPackage::SplitMultipoles(
97 const StaticSite& aps) const {
98
99 std::vector<QMPackage::MinimalMMCharge> multipoles_split;
100 // Calculate virtual charge positions
101 double a = options_.get("dipole_spacing").as<double>(); // this is in a0
102 double mag_d = aps.getDipole().norm(); // this is in e * a0
103 const Eigen::Vector3d dir_d = aps.getDipole().normalized();
104 const Eigen::Vector3d A = aps.getPos() + 0.5 * a * dir_d;
105 const Eigen::Vector3d B = aps.getPos() - 0.5 * a * dir_d;
106 double qA = mag_d / a;
107 double qB = -qA;
108 if (std::abs(qA) > 1e-12) {
109 multipoles_split.push_back(MinimalMMCharge(A, qA));
110 multipoles_split.push_back(MinimalMMCharge(B, qB));
111 }
112
113 if (aps.getRank() > 1) {
114 const Eigen::Matrix3d components = aps.CalculateCartesianMultipole();
115 Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> es;
116 es.computeDirect(components);
117 double a2 = 2 * a;
118 for (Index i = 0; i < 3; i++) {
119 double q = es.eigenvalues()[i] / (a2 * a2);
120 const Eigen::Vector3d vec1 =
121 aps.getPos() + 0.5 * a2 * es.eigenvectors().col(i);
122 const Eigen::Vector3d vec2 =
123 aps.getPos() - 0.5 * a2 * es.eigenvectors().col(i);
124 multipoles_split.push_back(MinimalMMCharge(vec1, q));
125 multipoles_split.push_back(MinimalMMCharge(vec2, q));
126 }
127 }
128 return multipoles_split;
129 }
130
GetLineAndSplit(std::ifstream & input_file,const std::string separators) const131 std::vector<std::string> QMPackage::GetLineAndSplit(
132 std::ifstream& input_file, const std::string separators) const {
133 std::string line;
134 tools::getline(input_file, line);
135 boost::trim(line);
136 return tools::Tokenizer(line, separators).ToVector();
137 }
138
139 } // namespace xtp
140 } // namespace votca
141