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> 39 void ClassicalSegment<T>::LoadFromFile(std::string filename) { Identify()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 } ParseOptions(const tools::Property & options)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 Evaluate(Topology & top)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); AddSteptoFilename(const std::string & filename,Index step)139 this->atomlist_.back().setMultipole(multipoles, rank); 140 readinmultipoles = 0; 141 } 142 } 143 } 144 this->calcPos(); 145 } 146 147 template <class T> 148 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> 157 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> 169 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> 190 std::string ClassicalSegment<T>::identify() const { 191 return ""; 192 } 193 194 template <> 195 std::string ClassicalSegment<PolarSite>::identify() const { 196 return "PolarSegment"; 197 } 198 template <> 199 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