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>();
      }
      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 meaningful local frame/rotation exists in this case (that is
    // exactly why this method, rather than MapMapAtomonMD, is being
    // used at all) -- the identity is the natural, correct transform
    // here, matching this method's own, existing, direct position
    // copy (a->getPos(), no rotation applied to it either).
    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);
    // Direction is translation-invariant by construction (unlike an
    // atom's own position), so only the rotation part of this
    // fragment's own MD->QM-template transform applies here -- no
    // shift_map2md/md_com involved at all, matching the same
    // "rotate the difference only" principle Atom::Rotate/
    // QMAtom::Rotate's own, existing implementation already uses for
    // ordinary positions.
    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);
  return map(seg, coordsfilename);
}
 
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);
 
  // Real, direct fix -- worked through directly with the user, given
  // this check's own earlier, exact equality requirement was
  // genuinely too strict: only the number of atoms actually
  // GENUINELY used later in this same function
  // (Result[id.first], further below, for every real qmatoms index)
  // needs to exist at all -- confirmed directly, by tracing every
  // real use of Result throughout this same function, that
  // correctness only ever genuinely requires each such index to be
  // in range, not that every single atom in the coordinate file gets
  // referenced by some qmatoms entry at all. A coordinate file with
  // MORE atoms than are actually referenced is a real, legitimate use
  // case (e.g. deliberately reusing the same, single template file
  // across several, differently-truncated mappings of the same
  // underlying fragment -- a real, direct case this exact check
  // originally, wrongly rejected) -- only genuinely too FEW atoms in
  // the file is a real, direct problem (a genuine, real
  // out-of-bounds risk below).
  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 real, actual safety guarantee needed -- every individual,
  // real qmatoms index referenced by the mapping file must itself be
  // in range, regardless of the file's own total atom count (the
  // check above alone does not guarantee this on its own, e.g. if the
  // mapping file's own qmatoms indices were not simply 0..N-1 at all
  // for some reason).
  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).");
    }
  }
 
  // Built incrementally, below, across all of this segment's own
  // fragments -- deliberately spans the WHOLE segment, not just one
  // fragment at a time, since a bonded partner can genuinely be in a
  // different fragment of the same segment (used only by
  // TransferBondedPartners, in a second pass after the main loop
  // below, once this is complete for the whole segment).
  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: now that md_id_to_map_atom is complete for the whole
  // segment, translate each atom's own, raw, MD-level bonded-partner
  // IDs into the corresponding QM-level ones. Deliberately a separate
  // pass, after the main loop above, rather than done inline within
  // it -- a bonded partner found while processing an earlier fragment
  // might belong to a later fragment not yet added to the lookup at
  // that point.
  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);
      }
    }
  }
 
  // Real, direct fix for a real, direct mistake in the earlier
  // relaxation of this function's own atom-count check (further
  // above): that fix correctly confirmed every real qmatoms index is
  // in range for Result[id.first] -- but missed that this function
  // used to simply `return Result;` directly, the WHOLE loaded
  // template file, not just the atoms actually referenced by
  // seginfo.mapatoms. Before that check's own earlier, strict
  // equality requirement was relaxed, this was harmless (equal
  // counts implicitly meant "every atom is referenced" too) -- but
  // once a coordinate file is legitimately allowed to have MORE
  // atoms than are referenced (the whole point of that fix), simply
  // returning Result directly leaks every genuinely unreferenced
  // template atom straight into the final, real QMMolecule/
  // StaticSegment/PolarSegment -- confirmed directly, this exact way,
  // from the user's own real, direct mapchecker output: an extra,
  // spurious H atom (single_thio.xyz's own 8th, deliberately
  // unreferenced atom) showing up in the actual mapped
  // "P3HTbackbone_1" segment, at exactly its own position.
  //
  // Fixed by building a real, separate, filtered AtomContainer
  // containing only the atoms genuinely referenced by
  // seginfo.mapatoms, in the SAME order the mapping file's own
  // qmatoms entries actually appear in (matching mdatoms' own,
  // already-established, canonical parallel order) -- but, critically
  // (a real, direct, second mistake caught and fixed here, worked
  // through directly with the user, in this same, real return-Result
  // bug): NOT simply copying each atom's own, original id_ forward.
  // The mapping file's own qmatoms entries are never required to be
  // sequential OR contiguous at all -- reordering a fully-referenced
  // set (e.g. 0:C 2:H 1:H 3:H, referencing every atom, just not in
  // file order) is explicitly a real, intended, legitimate use of
  // this same mapping mechanism, not merely the newer, adjacent
  // "leave some atoms out" case. Simply copying id_ forward would
  // silently give filtered[1] the atom whose own, original id was 2,
  // not 1 -- breaking the getId()-matches-position invariant this
  // whole codebase relies on elsewhere (e.g.
  // IPodCoupling::EvalJob's own fragment-A/B classification, which
  // compares atom.getId() directly against atom counts). Each copied
  // atom's own id is therefore explicitly renumbered, via setID(), to
  // match its own, real, new position within filtered -- not
  // reconstructed from scratch, since that would be unsafe,
  // generically, for StaticSite/PolarSite specifically (their own
  // real ClassicalSegment::LoadFromFile, classicalsegment.cc,
  // genuinely carries real multipole/rank data beyond plain
  // element/position, which a from-scratch reconstruction would
  // silently drop) -- setID() (added directly to StaticSite,
  // staticsite.h, matching QMAtom's own, already-existing setID()
  // naming exactly) preserves everything else about the atom exactly,
  // mutating only its own id.
  AtomContainer filtered(seg.getType(), seg.getId());
  for (const atom_id& id : seginfo.mapatoms) {
    mapAtom atom_copy = Result[id.first];
    atom_copy.setID(filtered.size());
    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