classicalsegment.cc

Go to the documentation of this file.

/*
 *            Copyright 2009-2020 The VOTCA Development Team
 *                       (http://www.votca.org)
 *
 *      Licensed under the Apache License, Version 2.0 (the "License")
 *
 * You may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *              http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */
 
// Standard includes
#include <algorithm>
#include <type_traits>
 
// Third party includes
#include <boost/format.hpp>
 
// VOTCA includes
#include <votca/tools/elements.h>
#include <votca/tools/getline.h>
#include <votca/tools/tokenizer.h>
 
// Local VOTCA includes
#include "votca/xtp/atomcontainer.h"
#include "votca/xtp/classicalsegment.h"
 
namespace votca {
namespace xtp {
 
// MPS files have a weird format positions can be in bohr or angstroem,
// multipoles are in q*bohr^k, with k rank of multipole and polarisabilities are
// in angstroem^3
template <class T>
void ClassicalSegment<T>::LoadFromFile(std::string filename) {
 
  std::string line;
  std::ifstream intt;
  intt.open(filename);
  double unit_conversion = tools::conv::ang2bohr;
 
  Index readinmultipoles = 0;
  Index numberofmultipoles = 0;
  Vector9d multipoles = Vector9d::Zero();
  Index rank = 0;
 
  if (!intt.is_open()) {
    throw std::runtime_error("File:" + filename + " could not be opened");
  }
  while (intt.good()) {
 
    tools::getline(intt, line);
    tools::Tokenizer toker(line, " \t");
    std::vector<std::string> split = toker.ToVector();
 
    if (!split.size() || split[0] == "!" || split[0].substr(0, 1) == "!") {
      continue;
    }
 
    // ! Interesting information here, e.g.
    // ! DCV2T opt
    // ! SP        RB3LYP          6-311+G(d,p)
    // Units bohr
    //
    // C          -4.2414603400   -3.8124751600    0.0017575736    Rank  2
    //  -0.3853409355
    //  -0.0002321905   0.2401559510   0.6602334308
    //  -0.7220625314   0.0004894995  -0.0003833545   0.4526409813  -0.50937399
    //  P 1.75
 
    // Units used
    if (split[0] == "Units") {
      std::string units = split[1];
      if (units != "bohr" && units != "angstrom") {
        throw std::runtime_error("Unit " + units + " in file " + filename +
                                 " not supported.");
      }
      if (units == "bohr") {
        unit_conversion = 1.0;
      }
    } else if (split.size() == 6) {
      // element,  position,  rank limit convert to bohr
      std::string name = split[0];
      Eigen::Vector3d pos;
      Index id = Index(this->atomlist_.size());
      pos[0] = boost::lexical_cast<double>(split[1]);
      pos[1] = boost::lexical_cast<double>(split[2]);
      pos[2] = boost::lexical_cast<double>(split[3]);
      rank = boost::lexical_cast<Index>(split[5]);
      numberofmultipoles = (rank + 1) * (rank + 1);
      multipoles = Vector9d::Zero();
      pos *= unit_conversion;
      this->atomlist_.push_back(T(id, name, pos));
    }
    // 'P', dipole polarizability
    else if (split[0] == "P") {
      Eigen::Matrix3d p1;
      if (split.size() == 7) {
        double pxx = boost::lexical_cast<double>(split[1]);
        double pxy = boost::lexical_cast<double>(split[2]);
        double pxz = boost::lexical_cast<double>(split[3]);
        double pyy = boost::lexical_cast<double>(split[4]);
        double pyz = boost::lexical_cast<double>(split[5]);
        double pzz = boost::lexical_cast<double>(split[6]);
        p1 << pxx, pxy, pxz, pxy, pyy, pyz, pxz, pyz, pzz;
      } else if (split.size() == 2) {
        double pxx = boost::lexical_cast<double>(split[1]);
        p1 = pxx * Eigen::Matrix3d::Identity();
      } else {
        throw std::runtime_error("Invalid line in " + filename + ": " + line);
      }
      double unit_conversion_3 = std::pow(tools::conv::ang2bohr, 3);
      p1 = p1 * unit_conversion_3;
      this->atomlist_.back().setpolarization(p1);
    }
    // Multipole lines
    else {
      // stay in bohr
      for (const std::string& entry : split) {
        double qXYZ = boost::lexical_cast<double>(entry);
        if (multipoles.size() < readinmultipoles) {
          throw std::runtime_error("ReadMpsFile: File" + filename +
                                   "is not properly formatted");
        }
        multipoles(readinmultipoles) = qXYZ;
        readinmultipoles++;
      }
      if (readinmultipoles == numberofmultipoles) {
        Eigen::Vector3d temp_dipoles = multipoles.segment<3>(1);
        // mps format for dipoles is z x y
        // we need x y z
        multipoles(1) = temp_dipoles(1);
        multipoles(2) = temp_dipoles(2);
        multipoles(3) = temp_dipoles(0);
        this->atomlist_.back().setMultipole(multipoles, rank);
        readinmultipoles = 0;
      }
    }
  }
  this->calcPos();
}
 
template <class T>
double ClassicalSegment<T>::CalcTotalQ() const {
  double Q = 0;
  for (const T& site : this->atomlist_) {
    Q += site.getCharge();
  }
  return Q;
}
 
template <class T>
Eigen::Vector3d ClassicalSegment<T>::CalcDipole() const {
  Eigen::Vector3d dipole = Eigen::Vector3d::Zero();
 
  Eigen::Vector3d CoM = this->getPos();
  for (const T& site : this->atomlist_) {
    dipole += (site.getPos() - CoM) * site.getCharge();
    dipole += site.getDipole();
  }
  return dipole;
}
 
template <class T>
Eigen::Vector3d ClassicalSegment<T>::CalcStaticDipole() const {
  Eigen::Vector3d dipole = Eigen::Vector3d::Zero();
 
  Eigen::Vector3d CoM = this->getPos();
  for (const T& site : this->atomlist_) {
    dipole += (site.getPos() - CoM) * site.getCharge();
    dipole += site.getStaticDipole();
  }
  return dipole;
}
 
template <class T>
void ClassicalSegment<T>::WriteMPS(std::string filename,
                                   std::string header) const {
 
  std::ofstream ofs;
  ofs.open(filename, std::ofstream::out);
  if (!ofs.is_open()) {
    throw std::runtime_error("Bad file handle: " + filename);
  }
 
  ofs << (boost::format("! GENERATED BY VOTCA::XTP::%1$s\n") % header);
  ofs << (boost::format("! N=%2$d Q[e]=%1$+1.7f\n") % CalcTotalQ() %
          this->size());
  ofs << boost::format("Units angstrom\n");
 
  for (const T& site : this->atomlist_) {
    ofs << site.WriteMpsLine("angstrom");
  }
  ofs.close();
}
 
template <class T>
void ClassicalSegment<T>::setChargeResponse(const ChargeResponse& cr) {
  if (cr.size() != this->size()) {
    throw std::runtime_error(
        "Segment " + this->getType() + " " + std::to_string(this->getId()) +
        " has " + std::to_string(this->size()) +
        " sites, its charge-response matrix " + std::to_string(cr.size()));
  }
  for (Index i = 0; i < this->size(); ++i) {
    if (cr.Elements()[std::size_t(i)] != this->atomlist_[i].getElement()) {
      throw std::runtime_error(
          "Segment " + this->getType() + " " + std::to_string(this->getId()) +
          ": site " + std::to_string(i) + " is " +
          this->atomlist_[i].getElement() +
          ", but the charge-response matrix has " +
          cr.Elements()[std::size_t(i)] + " there (atom order must agree)");
    }
  }
  chargeresponse_ = cr;
}
 
template <class T>
void ClassicalSegment<T>::WriteToCpt(CheckpointWriter& w) const {
  AtomContainer<T>::WriteToCpt(w);
  if (hasChargeResponse()) {
    CheckpointWriter v = w.openChild("chargeresponse");
    chargeresponse_.WriteToCpt(v);
    if constexpr (std::is_same_v<T, PolarSite>) {
      Eigen::VectorXd dq(this->size());
      for (Index i = 0; i < this->size(); ++i) {
        dq(i) = this->atomlist_[i].Induced_Charge();
      }
      v(dq, "induced_charges");
    }
  }
}
 
template <class T>
void ClassicalSegment<T>::ReadFromCpt(CheckpointReader& r) {
  AtomContainer<T>::ReadFromCpt(r);
  ReadChargeResponseFromCpt(r);
}
 
template <class T>
void ClassicalSegment<T>::ReadChargeResponseFromCpt(CheckpointReader& r) {
  chargeresponse_ = ChargeResponse();
  // optional: checkpoints without it (and all older ones) are Thole segments
  const std::vector<std::string> children = r.getChildGroupNames();
  if (std::find(children.begin(), children.end(), "chargeresponse") !=
      children.end()) {
    CheckpointReader v = r.openChild("chargeresponse");
    ChargeResponse cr;
    cr.ReadFromCpt(v);
    setChargeResponse(cr);
    if constexpr (std::is_same_v<T, PolarSite>) {
      Eigen::VectorXd dq;
      v(dq, "induced_charges");
      if (dq.size() != this->size()) {
        throw std::runtime_error("Segment " + this->getType() +
                                 ": induced charges in the checkpoint do not "
                                 "match its sites");
      }
      for (Index i = 0; i < this->size(); ++i) {
        this->atomlist_[i].setInduced_Charge(dq(i));
      }
    }
  }
}
 
template <class T>
std::string ClassicalSegment<T>::identify() const {
  return "";
}
 
template <>
std::string ClassicalSegment<PolarSite>::identify() const {
  return "PolarSegment";
}
template <>
std::string ClassicalSegment<StaticSite>::identify() const {
  return "StaticSegment";
}
 
// additional access functions and methods for Ewald
template <class T>
bool ClassicalSegment<T>::IsCharged() const {
  for (const T& site : this->atomlist_) {
    if (std::abs(site.getCharge()) > 1e-4) {
      return true;
    }
  }
  return false;
}
 
template <class T>
bool ClassicalSegment<T>::IsPolarizable() const {
  return false;
}
 
template <>
bool ClassicalSegment<PolarSite>::IsPolarizable() const {
  for (const PolarSite& site : this->atomlist_) {
 
    (site.getpolarization().array().abs() > 1e-4).any();
    if ((site.getpolarization().array().abs() > 1e-4).any()) {
      return true;
    }
  }
  return false;
}
 
template <>
bool ClassicalSegment<StaticSite>::IsPolarizable() const {
  return false;
}
 
template class ClassicalSegment<PolarSite>;
template class ClassicalSegment<StaticSite>;
 
}  // namespace xtp
}  // namespace votca