segmentmapper.cc

Go to the documentation of this file.

/*
 *            Copyright 2016 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 <regex>
 
// Local VOTCA includes
#include "votca/xtp/segmentmapper.h"
 
namespace votca {
namespace xtp {
 
template <class AtomContainer>
SegmentMapper<AtomContainer>::SegmentMapper(Logger& log) : log_(log) {
  FillMap();
}
 
template <class AtomContainer>
AtomContainer SegmentMapper<AtomContainer>::map(const Segment& seg,
                                                const SegId& segid) const {
  if (segid.hasFile()) {
    std::string filename = segid.FileName();
    filename = std::regex_replace(filename, std::regex("\\$SEGID"),
                                  std::to_string(seg.getId()));
    filename =
        std::regex_replace(filename, std::regex("\\$SEGNAME"), seg.getType());
    return map(seg, filename);
  } else {
    QMState state = segid.getQMState();
    return map(seg, state);
  }
}
 
template <class AtomContainer>
template <typename T>
Eigen::Vector3d SegmentMapper<AtomContainer>::CalcWeightedPos(
    const std::vector<double>& weights, const T& atoms) const {
  Eigen::Vector3d map_pos = Eigen::Vector3d::Zero();
  double map_weight = 0.0;
  for (Index i = 0; i < Index(atoms.size()); i++) {
    map_pos += atoms[i]->getPos() * weights[i];
    map_weight += weights[i];
  }
  return map_pos = map_pos / map_weight;
}
 
template <class AtomContainer>
std::vector<double> SegmentMapper<AtomContainer>::getWeights(
    const tools::Property& frag) const {
 
  std::vector<double> weights;
  if (frag.exists(mapatom_xml_.at("weights"))) {
    weights =
        frag.get(mapatom_xml_.at("weights")).template as<std::vector<double>>();
  } else if (frag.exists("weights")) {
    weights = frag.get("weights").template as<std::vector<double>>();
  } else {
    XTP_LOG(Log::error, log_) << " Did not find weights for fragment "
                              << frag.get("name").as<std::string>()
                              << " Using atomic masses" << std::flush;
    std::string frags =
        frag.get(mapatom_xml_.at("atoms")).template as<std::string>();
    tools::Tokenizer tok_atoms(frags, " \t\n");
    std::vector<std::string> atom_strings = tok_atoms.ToVector();
    tools::Elements e;
    for (auto a_string : atom_strings) {
      tools::Tokenizer tok_atom(a_string, ":");
      std::vector<std::string> entries = tok_atom.ToVector();
      if (entries.size() != 2) {
        throw std::runtime_error("Cannot get weight from Element for " +
                                 a_string);
      }
 
      double weight = e.getMass(entries[1]);
      XTP_LOG(Log::info, log_) << entries[1] << ":" << weight << " ";
      weights.push_back(weight);
    }
    XTP_LOG(Log::error, log_) << std::endl;
  }
 
  return weights;
}
 
template <class AtomContainer>
void SegmentMapper<AtomContainer>::ParseFragment(Seginfo& seginfo,
                                                 const tools::Property& frag) {
  std::vector<std::string> map_atoms =
      frag.get(mapatom_xml_["atoms"]).template as<std::vector<std::string>>();
  std::vector<std::string> md_atoms =
      frag.get("mdatoms").as<std::vector<std::string>>();
 
  if (md_atoms.size() != map_atoms.size()) {
    throw std::runtime_error(
        "Mapping for segment " + seginfo.segname + " fragment " +
        frag.get("name").as<std::string>() +
        " does not have same numbers of md and " + mapatom_xml_["atoms"] +
        "."
        "If you want to leave a qmatom out, place a ':' instead");
  }
 
  std::vector<double> weights = getWeights(frag);
 
  if (md_atoms.size() != weights.size()) {
    throw std::runtime_error("Mapping for segment " + seginfo.segname +
                             " fragment " + frag.get("name").as<std::string>() +
                             " does not have same numbers of md and weights. "
                             "If you want to leave a " +
                             mapatom_xml_["name"] +
                             " out, place a '0' instead");
  }
 
  FragInfo mapfragment;
  std::vector<Index> mapatom_ids;
 
  for (Index i = 0; i < Index(map_atoms.size()); i++) {
    const std::string& map_string = map_atoms[i];
    const std::string& md_string = md_atoms[i];
    const double& weight = weights[i];
    atom_id md_result = StringToMDIndex(md_string);
    seginfo.mdatoms.push_back(md_result.first);
 
    if (map_string == ":") {
      continue;
    }
    atom_id map_result = StringToMapIndex(map_string);
    mapatom_ids.push_back(map_result.first);
    seginfo.mapatoms.push_back(map_result);
    if (Atom::GetElementFromString(md_result.second) != map_result.second) {
      XTP_LOG(Log::error, log_)
          << "WARNING: mdatom'" << md_result.second << "' and "
          << mapatom_xml_["name"] << " '" << map_result.second
          << "' do not have same element" << std::flush;
    }
    mapfragment.mapatom_ids.push_back(map_result);
    mapfragment.weights.push_back(weight);
    mapfragment.mdatom_ids.push_back(md_result);
  }
 
  std::vector<Index> frame =
      tools::Tokenizer(getFrame(frag), " \t\n").ToVector<Index>();
  if (frame.size() > 3) {
    throw std::runtime_error(
        "Local frame for segment " + seginfo.segname + " fragment " +
        frag.get("name").as<std::string>() +
        " has more than 3 atoms, please specify only up to three atoms");
  } else if (frame.empty()) {
    throw std::runtime_error("No local frame for segment " + seginfo.segname +
                             " fragment " + frag.get("name").as<std::string>() +
                             " specified");
  }
  for (Index atomid : frame) {
    if (std::find(mapatom_ids.begin(), mapatom_ids.end(), atomid) ==
        mapatom_ids.end()) {
      throw std::runtime_error("Atom " + std::to_string(atomid) +
                               " in local frame cannot be found in " +
                               mapatom_xml_["atoms"] + ".");
    }
  }
 
  mapfragment.map_local_frame = frame;
  seginfo.fragments.push_back(mapfragment);
}
 
template <class AtomContainer>
void SegmentMapper<AtomContainer>::LoadMappingFile(const std::string& mapfile) {
  tools::Property topology_map;
  topology_map.LoadFromXML(mapfile);
 
  std::string molkey = "topology.molecules.molecule";
  std::vector<tools::Property*> molecules = topology_map.Select(molkey);
  std::string segkey = "segments.segment";
  for (tools::Property* mol : molecules) {
    std::vector<tools::Property*> segments = mol->Select(segkey);
    for (tools::Property* seg : segments) {
      Seginfo seginfo;
 
      std::string coordfile_key = mapatom_xml_["coords"] + "_*";
      std::vector<tools::Property*> files = seg->Select(coordfile_key);
      for (tools::Property* file : files) {
        seginfo.coordfiles[file->name()] = file->as<std::string>();
      }
      for (tools::Property* file : seg->Select("chargeresponse_*")) {
        seginfo.crfiles[file->name()] = file->as<std::string>();
      }
      std::string segname = seg->get("name").as<std::string>();
      seginfo.segname = segname;
      std::string fragkey = "fragments.fragment";
 
      seginfo.map2md = seg->get("map2md").as<bool>();
 
      std::vector<tools::Property*> fragments = seg->Select(fragkey);
      for (tools::Property* frag : fragments) {
        ParseFragment(seginfo, *frag);
      }
 
      Index map_atom_min_id =
          std::min_element(seginfo.mapatoms.begin(), seginfo.mapatoms.end(),
                           [](const atom_id& a, const atom_id& b) {
                             return a.first < b.first;
                           })
              ->first;
      if (map_atom_min_id != 0) {
        throw std::runtime_error(
            mapatom_xml_["atoms"] + " for segment " + seginfo.segname +
            " do not start at zero index. Each segment "
            "should have its own coordinate file. If you use an old ctp "
            "mapping file run 'xtp_update_mapfile' on it.");
      }
 
      seginfo.minmax = CalcAtomIdRange(seginfo.mdatoms);
      segment_info_[segname] = seginfo;
    }
  }
}
 
template <class AtomContainer>
std::pair<Index, std::string> SegmentMapper<AtomContainer>::StringToMapIndex(
    const std::string& map_string) const {
  tools::Tokenizer tok(map_string, ":");
  std::vector<std::string> result = tok.ToVector();
  if (result.size() != 2) {
    throw std::runtime_error("Entry " + map_string +
                             " is not properly formatted.");
  }
  return std::pair<Index, std::string>(std::stoi(result[0]), result[1]);
}
template <class AtomContainer>
std::pair<Index, std::string> SegmentMapper<AtomContainer>::StringToMDIndex(
    const std::string& md_string) const {
  tools::Tokenizer tok(md_string, ":");
  std::vector<std::string> result = tok.ToVector();
  if (result.size() != 3) {
    throw std::runtime_error("Entry " + md_string +
                             " is not properly formatted.");
  }
  Index atomid = 0;
  try {
    atomid = std::stoi(result[2]);
  } catch (std::invalid_argument&) {
    throw std::runtime_error("Atom entry " + md_string +
                             " is not well formatted");
  }
  return std::pair<Index, std::string>(atomid, result[1]);
}
 
template <class AtomContainer>
std::pair<Index, Index> SegmentMapper<AtomContainer>::CalcAtomIdRange(
    const std::vector<Index>& seg) const {
  Index max_res_id = *std::max_element(seg.begin(), seg.end());
  Index min_res_id = *std::min_element(seg.begin(), seg.end());
  return std::pair<Index, Index>(min_res_id, max_res_id);
}
template <class AtomContainer>
std::pair<Index, Index> SegmentMapper<AtomContainer>::CalcAtomIdRange(
    const Segment& seg) const {
  Index max_res_id = std::max_element(seg.begin(), seg.end(),
                                      [](const Atom& a, const Atom& b) {
                                        return a.getId() < b.getId();
                                      })
                         ->getId();
 
  Index min_res_id = std::min_element(seg.begin(), seg.end(),
                                      [](const Atom& a, const Atom& b) {
                                        return a.getId() < b.getId();
                                      })
                         ->getId();
  return std::pair<Index, Index>(min_res_id, max_res_id);
}
 
template <class AtomContainer>
void SegmentMapper<AtomContainer>::PlaceMapAtomonMD(
    const std::vector<mapAtom*>& fragment_mapatoms,
    const std::vector<const Atom*>& fragment_mdatoms) const {
  for (Index i = 0; i < Index(fragment_mapatoms.size()); i++) {
    const Atom* a = fragment_mdatoms[i];
    mapAtom* b = fragment_mapatoms[i];
    b->setPos(a->getPos());
    // No local frame is used here, positions are copied unrotated, so the
    // bond direction is transferred with the identity.
    TransferExternalBondDirection(b, a, Eigen::Matrix3d::Identity());
  }
}
 
template <class AtomContainer>
Index SegmentMapper<AtomContainer>::FindVectorIndexFromAtomId(
    Index atomid, const std::vector<mapAtom*>& fragment_mapatoms) const {
  Index i = 0;
  for (; i < Index(fragment_mapatoms.size()); i++) {
    if (fragment_mapatoms[i]->getId() == atomid) {
      break;
    }
  }
  return i;
}
template <class AtomContainer>
void SegmentMapper<AtomContainer>::MapMapAtomonMD(
    const FragInfo& frag, const std::vector<mapAtom*>& fragment_mapatoms,
    const std::vector<const Atom*>& fragment_mdatoms) const {
  std::vector<Eigen::Vector3d> local_map_frame;
  std::vector<Eigen::Vector3d> local_md_frame;
  for (Index id : frag.map_local_frame) {
    Index i = FindVectorIndexFromAtomId(id, fragment_mapatoms);
    local_map_frame.push_back(fragment_mapatoms[i]->getPos());
    local_md_frame.push_back(fragment_mdatoms[i]->getPos());
  }
 
  Index symmetry = frag.map_local_frame.size();
  Eigen::Vector3d map_com = CalcWeightedPos(frag.weights, fragment_mapatoms);
  Eigen::Vector3d md_com = CalcWeightedPos(frag.weights, fragment_mdatoms);
 
  Eigen::Vector3d shift_map2md = md_com - map_com;
 
  Eigen::Matrix3d rot_map = Eigen::Matrix3d::Identity();
  Eigen::Matrix3d rot_md = Eigen::Matrix3d::Identity();
 
  // building local frame
  if (symmetry > 1) {
    // middle atom is the origin
    Eigen::Vector3d x_map = local_map_frame[0] - local_map_frame[1];
    Eigen::Vector3d x_md = local_md_frame[0] - local_md_frame[1];
    Eigen::Vector3d y_map;
    Eigen::Vector3d y_md;
    Eigen::Vector3d z_map;
    Eigen::Vector3d z_md;
    if (symmetry == 3) {
      y_map = local_map_frame[2] - local_map_frame[1];
      y_md = local_md_frame[2] - local_md_frame[1];
      z_map = x_map.cross(y_map);
      z_md = x_md.cross(y_md);
      y_map = z_map.cross(x_map);
      y_md = z_md.cross(x_md);
 
    } else {
      Eigen::Vector3d unit = Eigen::Vector3d::UnitX();
      if (std::abs(unit.dot(x_map) / x_map.norm()) < 1e-6) {
        unit = Eigen::Vector3d::UnitY();
      }
      y_map = (x_map.cross(unit));
      if (std::abs(unit.dot(x_md) / x_md.norm()) < 1e-6) {
        unit = Eigen::Vector3d::UnitX();
      }
      y_md = (x_md.cross(unit));
      z_map = x_map.cross(y_map);
      z_md = x_md.cross(y_md);
    }
 
    if (x_map.squaredNorm() < 1e-18 || y_map.squaredNorm() < 1e-18 ||
        z_map.squaredNorm() < 1e-18) {
      throw std::runtime_error(
          mapatom_xml_.at("tag") +
          " basis vectors are very small, choose different local basis");
    }
    rot_map.col(0) = x_map.normalized();
    rot_map.col(1) = y_map.normalized();
    rot_map.col(2) = z_map.normalized();
 
    if (x_md.squaredNorm() < 1e-18 || y_md.squaredNorm() < 1e-18 ||
        z_md.squaredNorm() < 1e-18) {
      throw std::runtime_error(
          "MD basis vectors are very small, choose different local basis");
    }
    rot_md.col(0) = x_md.normalized();
    rot_md.col(1) = y_md.normalized();
    rot_md.col(2) = z_md.normalized();
  }
  Eigen::Matrix3d rotateMAP2MD = rot_md * rot_map.transpose();
  for (Index i = 0; i < Index(fragment_mapatoms.size()); i++) {
    mapAtom* atom = fragment_mapatoms[i];
    if (getRank(*atom) > 0 && symmetry < 3) {
      throw std::runtime_error(
          "Local frame has less than 3 atoms, thus higher rank multipoles "
          "cannot be mapped.");
    }
    atom->Translate(shift_map2md);
    atom->Rotate(rotateMAP2MD, md_com);
    // A direction is translation invariant, so only the rotation applies.
    TransferExternalBondDirection(atom, fragment_mdatoms[i], rotateMAP2MD);
  }
}
template <class AtomContainer>
AtomContainer SegmentMapper<AtomContainer>::map(const Segment& seg,
                                                QMState state) const {
  if (segment_info_.count(seg.getType()) == 0) {
    throw std::runtime_error(
        "Could not find a Segment of name: " + seg.getType() + " in mapfile.");
  }
  Seginfo seginfo = segment_info_.at(seg.getType());
  std::string coordsfiletag =
      mapatom_xml_.at("coords") + "_" + state.Type().ToString();
  if (seginfo.coordfiles.count(coordsfiletag) == 0) {
    throw std::runtime_error("Could not find a coordinate file for " +
                             seg.getType() +
                             " id:" + std::to_string(seg.getId()) +
                             " segment/state: " + coordsfiletag);
  }
  std::string coordsfilename = seginfo.coordfiles.at(coordsfiletag);
  AtomContainer result = map(seg, coordsfilename);
  AttachChargeResponse(result, seginfo, state);
  return result;
}
 
template <>
void SegmentMapper<PolarSegment>::AttachChargeResponse(
    PolarSegment& result, const Seginfo& seginfo, const QMState& state) const {
  const std::string tag = "chargeresponse_" + state.Type().ToString();
  if (seginfo.crfiles.count(tag) == 0) {
    return;
  }
  const std::string& filename = seginfo.crfiles.at(tag);
  ChargeResponse cr;
  cr.LoadFromFile(filename);
  // The mapped segment holds the template atoms seginfo.mapatoms in that
  // order; Pi must follow them. Every atom must be mapped: Pi of a part of
  // a molecule does not conserve charge.
  std::vector<Index> order;
  for (const atom_id& id : seginfo.mapatoms) {
    if (id.first >= cr.size() ||
        cr.Elements()[std::size_t(id.first)] != id.second) {
      throw std::runtime_error(
          "Segment " + seginfo.segname + ": " + mapatom_xml_.at("atoms") +
          " entry " + std::to_string(id.first) + ":" + id.second +
          " does not match the charge-response file " + filename +
          " (same atoms in the same order as the multipoles file needed)");
    }
    order.push_back(id.first);
  }
  if (Index(order.size()) != cr.size()) {
    throw std::runtime_error(
        "Segment " + seginfo.segname + " maps " + std::to_string(order.size()) +
        " of the " + std::to_string(cr.size()) +
        " atoms of its charge-response file " + filename +
        "; the charge-response model needs the whole molecule");
  }
  result.setChargeResponse(cr.Permuted(order));
  XTP_LOG(Log::info, log_) << "Segment " << seginfo.segname << " "
                           << result.getId() << ": charge response from "
                           << filename << std::flush;
}
 
template <class AtomContainer>
AtomContainer SegmentMapper<AtomContainer>::map(
    const Segment& seg, const std::string& coordfilename) const {
 
  if (segment_info_.count(seg.getType()) == 0) {
    throw std::runtime_error(
        "Could not find a Segment of name: " + seg.getType() + " in mapfile.");
  }
  Seginfo seginfo = segment_info_.at(seg.getType());
  if (Index(seginfo.mdatoms.size()) != seg.size()) {
    throw std::runtime_error(
        "Segment '" + seg.getType() +
        "' does not contain the same number of atoms as mapping file: " +
        std::to_string(seginfo.mdatoms.size()) + " vs. " +
        std::to_string(seg.size()));
  }
 
  std::pair<Index, Index> minmax_map = seginfo.minmax;
  std::pair<Index, Index> minmax = CalcAtomIdRange(seg);
 
  if ((minmax_map.first - minmax_map.second) !=
      (minmax.first - minmax.second)) {
    throw std::runtime_error("AtomId range for segment " + seg.getType() + ":" +
                             std::to_string(seg.getId()) +
                             " and the mapping do not agree: Segment[" +
                             std::to_string(minmax.first) + "," +
                             std::to_string(minmax.second) + "] Map[" +
                             std::to_string(minmax_map.first) + "," +
                             std::to_string(minmax_map.second) + "]");
  }
 
  Index atomidoffset = minmax.first - minmax_map.first;
 
  AtomContainer Result(seg.getType(), seg.getId());
  Result.LoadFromFile(coordfilename);
 
  // The coordinate file may contain more atoms than the mapping references
  // (e.g. one template file reused for differently truncated mappings);
  // only too few atoms is an error.
  if (Result.size() < Index(seginfo.mapatoms.size())) {
    throw std::runtime_error(
        mapatom_xml_.at("tag") + "Segment '" + seg.getType() +
        "' does not contain enough atoms for the mapping file: needs at "
        "least " +
        std::to_string(seginfo.mapatoms.size()) + ", has " +
        std::to_string(Result.size()));
  }
  // The count check above does not guarantee that each referenced index is
  // in range, since the indices need not be 0..N-1.
  for (const atom_id& id : seginfo.mapatoms) {
    if (id.first < 0 || id.first >= Result.size()) {
      throw std::runtime_error(
          mapatom_xml_.at("tag") + "Segment '" + seg.getType() +
          "' references " + mapatom_xml_.at("atoms") + " index " +
          std::to_string(id.first) +
          ", which is out of range for the coordinate file '" + coordfilename +
          "' (" + std::to_string(Result.size()) + " atoms).");
    }
  }
 
  // MD atom ID -> mapped atom for the whole segment (a bonded partner may be
  // in another fragment); used by TransferBondedPartners below.
  std::map<Index, mapAtom*> md_id_to_map_atom;
 
  for (FragInfo& frag : seginfo.fragments) {
    for (atom_id& id : frag.mdatom_ids) {
      id.first += atomidoffset;
    }
 
    std::vector<mapAtom*> fragment_mapatoms;
    for (const atom_id& id : frag.mapatom_ids) {
      if (id.second != Result[id.first].getElement()) {
        throw std::runtime_error(
            "Element of mapping atom " + std::to_string(id.first) + ":" +
            id.second + " does not agree with Element of parsed Atom " +
            Result[id.first].getElement());
      }
      fragment_mapatoms.push_back(&Result[id.first]);
    }
    std::vector<const Atom*> fragment_mdatoms;
    for (const atom_id& id : frag.mdatom_ids) {
      const Atom* atom = seg.getAtom(id.first);
      if (atom == nullptr) {
        throw std::runtime_error(
            "Could not find an atom with name:" + id.second + "id" +
            std::to_string(id.first) + " in segment " + seg.getType());
      }
      fragment_mdatoms.push_back(atom);
    }
 
    for (Index i = 0; i < Index(fragment_mdatoms.size()); i++) {
      md_id_to_map_atom[fragment_mdatoms[i]->getId()] = fragment_mapatoms[i];
    }
 
    if (seginfo.map2md) {
      PlaceMapAtomonMD(fragment_mapatoms, fragment_mdatoms);
    } else {
      MapMapAtomonMD(frag, fragment_mapatoms, fragment_mdatoms);
    }
  }
 
  // Second pass, once md_id_to_map_atom is complete: translate MD-level
  // bonded-partner IDs to mapped IDs. Partners in later fragments are not in
  // the lookup during the first loop.
  for (FragInfo& frag : seginfo.fragments) {
    for (const atom_id& id : frag.mdatom_ids) {
      const Atom* md_atom = seg.getAtom(id.first);
      if (md_atom == nullptr) {
        continue;
      }
      auto it = md_id_to_map_atom.find(id.first);
      if (it != md_id_to_map_atom.end()) {
        TransferBondedPartners(it->second, md_atom, md_id_to_map_atom);
      }
    }
  }
 
  // Return only the atoms referenced by seginfo.mapatoms, in mapping-file
  // order (parallel to mdatoms); unreferenced template atoms must not leak
  // into the result. The qmatoms entries may be reordered (e.g.
  // 0:C 2:H 1:H 3:H), so each copy is renumbered with setID() to keep
  // getId() equal to the position in the container, which e.g. the
  // fragment A/B classification in IPodCoupling::EvalJob relies on.
  // Atoms are copied rather than rebuilt to keep multipole data of
  // StaticSite/PolarSite. Bonded-partner IDs refer to template atom IDs and
  // are rewritten to the new numbering.
  std::map<Index, Index> old_to_new;
  for (Index k = 0; k < Index(seginfo.mapatoms.size()); k++) {
    old_to_new[Result[seginfo.mapatoms[k].first].getId()] = k;
  }
  AtomContainer filtered(seg.getType(), seg.getId());
  for (const atom_id& id : seginfo.mapatoms) {
    mapAtom atom_copy = Result[id.first];
    atom_copy.setID(filtered.size());
    RemapBondedPartners(atom_copy, old_to_new);
    filtered.push_back(atom_copy);
  }
  filtered.calcPos();
  return filtered;
}
 
template class SegmentMapper<QMMolecule>;
template class SegmentMapper<StaticSegment>;
template class SegmentMapper<PolarSegment>;
 
}  // namespace xtp
}  // namespace votca