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/format.hpp>
22
23 // VOTCA includes
24 #include <votca/tools/elements.h>
25 #include <votca/tools/getline.h>
26 #include <votca/tools/tokenizer.h>
27
28 // Local VOTCA includes
29 #include "votca/xtp/atomcontainer.h"
30 #include "votca/xtp/classicalsegment.h"
31
32 namespace votca {
33 namespace xtp {
34
35 // MPS files have a weird format positions can be in bohr or angstroem,
36 // multipoles are in q*bohr^k, with k rank of multipole and polarisabilities are
37 // in angstroem^3
38 template <class T>
LoadFromFile(std::string filename)39 void ClassicalSegment<T>::LoadFromFile(std::string filename) {
40
41 std::string line;
42 std::ifstream intt;
43 intt.open(filename);
44 double unit_conversion = tools::conv::ang2bohr;
45
46 Index readinmultipoles = 0;
47 Index numberofmultipoles = 0;
48 Vector9d multipoles = Vector9d::Zero();
49 Index rank = 0;
50
51 if (!intt.is_open()) {
52 throw std::runtime_error("File:" + filename + " could not be opened");
53 }
54 while (intt.good()) {
55
56 tools::getline(intt, line);
57 tools::Tokenizer toker(line, " \t");
58 std::vector<std::string> split = toker.ToVector();
59
60 if (!split.size() || split[0] == "!" || split[0].substr(0, 1) == "!") {
61 continue;
62 }
63
64 // ! Interesting information here, e.g.
65 // ! DCV2T opt
66 // ! SP RB3LYP 6-311+G(d,p)
67 // Units bohr
68 //
69 // C -4.2414603400 -3.8124751600 0.0017575736 Rank 2
70 // -0.3853409355
71 // -0.0002321905 0.2401559510 0.6602334308
72 // -0.7220625314 0.0004894995 -0.0003833545 0.4526409813 -0.50937399
73 // P 1.75
74
75 // Units used
76 if (split[0] == "Units") {
77 std::string units = split[1];
78 if (units != "bohr" && units != "angstrom") {
79 throw std::runtime_error("Unit " + units + " in file " + filename +
80 " not supported.");
81 }
82 if (units == "bohr") {
83 unit_conversion = 1.0;
84 }
85 } else if (split.size() == 6) {
86 // element, position, rank limit convert to bohr
87 std::string name = split[0];
88 Eigen::Vector3d pos;
89 Index id = Index(this->atomlist_.size());
90 pos[0] = boost::lexical_cast<double>(split[1]);
91 pos[1] = boost::lexical_cast<double>(split[2]);
92 pos[2] = boost::lexical_cast<double>(split[3]);
93 rank = boost::lexical_cast<Index>(split[5]);
94 numberofmultipoles = (rank + 1) * (rank + 1);
95 multipoles = Vector9d::Zero();
96 pos *= unit_conversion;
97 this->atomlist_.push_back(T(id, name, pos));
98 }
99 // 'P', dipole polarizability
100 else if (split[0] == "P") {
101 Eigen::Matrix3d p1;
102 if (split.size() == 7) {
103 double pxx = boost::lexical_cast<double>(split[1]);
104 double pxy = boost::lexical_cast<double>(split[2]);
105 double pxz = boost::lexical_cast<double>(split[3]);
106 double pyy = boost::lexical_cast<double>(split[4]);
107 double pyz = boost::lexical_cast<double>(split[5]);
108 double pzz = boost::lexical_cast<double>(split[6]);
109 p1 << pxx, pxy, pxz, pxy, pyy, pyz, pxz, pyz, pzz;
110 } else if (split.size() == 2) {
111 double pxx = boost::lexical_cast<double>(split[1]);
112 p1 = pxx * Eigen::Matrix3d::Identity();
113 } else {
114 throw std::runtime_error("Invalid line in " + filename + ": " + line);
115 }
116 double unit_conversion_3 = std::pow(tools::conv::ang2bohr, 3);
117 p1 = p1 * unit_conversion_3;
118 this->atomlist_.back().setpolarization(p1);
119 }
120 // Multipole lines
121 else {
122 // stay in bohr
123 for (const std::string& entry : split) {
124 double qXYZ = boost::lexical_cast<double>(entry);
125 if (multipoles.size() < readinmultipoles) {
126 throw std::runtime_error("ReadMpsFile: File" + filename +
127 "is not properly formatted");
128 }
129 multipoles(readinmultipoles) = qXYZ;
130 readinmultipoles++;
131 }
132 if (readinmultipoles == numberofmultipoles) {
133 Eigen::Vector3d temp_dipoles = multipoles.segment<3>(1);
134 // mps format for dipoles is z x y
135 // we need x y z
136 multipoles(1) = temp_dipoles(1);
137 multipoles(2) = temp_dipoles(2);
138 multipoles(3) = temp_dipoles(0);
139 this->atomlist_.back().setMultipole(multipoles, rank);
140 readinmultipoles = 0;
141 }
142 }
143 }
144 this->calcPos();
145 }
146
147 template <class T>
CalcTotalQ() const148 double ClassicalSegment<T>::CalcTotalQ() const {
149 double Q = 0;
150 for (const T& site : this->atomlist_) {
151 Q += site.getCharge();
152 }
153 return Q;
154 }
155
156 template <class T>
CalcDipole() const157 Eigen::Vector3d ClassicalSegment<T>::CalcDipole() const {
158 Eigen::Vector3d dipole = Eigen::Vector3d::Zero();
159
160 Eigen::Vector3d CoM = this->getPos();
161 for (const T& site : this->atomlist_) {
162 dipole += (site.getPos() - CoM) * site.getCharge();
163 dipole += site.getDipole();
164 }
165 return dipole;
166 }
167
168 template <class T>
WriteMPS(std::string filename,std::string header) const169 void ClassicalSegment<T>::WriteMPS(std::string filename,
170 std::string header) const {
171
172 std::ofstream ofs;
173 ofs.open(filename, std::ofstream::out);
174 if (!ofs.is_open()) {
175 throw std::runtime_error("Bad file handle: " + filename);
176 }
177
178 ofs << (boost::format("! GENERATED BY VOTCA::XTP::%1$s\n") % header);
179 ofs << (boost::format("! N=%2$d Q[e]=%1$+1.7f\n") % CalcTotalQ() %
180 this->size());
181 ofs << boost::format("Units angstrom\n");
182
183 for (const T& site : this->atomlist_) {
184 ofs << site.WriteMpsLine("angstrom");
185 }
186 ofs.close();
187 }
188
189 template <class T>
identify() const190 std::string ClassicalSegment<T>::identify() const {
191 return "";
192 }
193
194 template <>
identify() const195 std::string ClassicalSegment<PolarSite>::identify() const {
196 return "PolarSegment";
197 }
198 template <>
identify() const199 std::string ClassicalSegment<StaticSite>::identify() const {
200 return "StaticSegment";
201 }
202
203 template class ClassicalSegment<PolarSite>;
204 template class ClassicalSegment<StaticSite>;
205
206 } // namespace xtp
207 } // namespace votca
208