segmentmapper.h

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.
 *
 */
 
#pragma once
#ifndef VOTCA_XTP_SEGMENTMAPPER_H
#define VOTCA_XTP_SEGMENTMAPPER_H
 
// Standard includes
#include <type_traits>
 
// VOTCA includes
#include <votca/csg/pdbwriter.h>
#include <votca/tools/property.h>
 
// Local VOTCA includes
#include "classicalsegment.h"
#include "logger.h"
#include "qmmolecule.h"
#include "segid.h"
#include "topology.h"
 
namespace votca {
namespace xtp {
template <class AtomContainer>
class SegmentMapper {
 public:
  SegmentMapper(Logger& log);
 
  void LoadMappingFile(const std::string& mapfile);
 
  AtomContainer map(const Segment& seg, const SegId& segid) const;
 
  AtomContainer map(const Segment& seg, QMState state) const;
 
  AtomContainer map(const Segment& seg, const std::string& coordfilename) const;
 
 private:
  using mapAtom = typename AtomContainer::Atom_Type;
 
  using atom_id = std::pair<Index, std::string>;
 
  struct FragInfo {
    std::vector<double> weights;
    std::vector<atom_id> mapatom_ids;
    std::vector<atom_id> mdatom_ids;
    std::vector<Index> map_local_frame;
  };
 
  struct Seginfo {
    std::pair<Index, Index> minmax;
    std::vector<Index> mdatoms;
    std::vector<FragInfo> fragments;
    bool map2md;
    std::string segname;
    std::vector<double> weights;
    std::vector<atom_id> mapatoms;
    std::map<std::string, std::string> coordfiles;
  };
  std::map<std::string, std::string> mapatom_xml_;
  std::map<std::string, Seginfo> segment_info_;
 
  Index FindVectorIndexFromAtomId(
      Index atomid, const std::vector<mapAtom*>& fragment_mapatoms) const;
 
  void ParseFragment(Seginfo& seginfo, const tools::Property& frag);
 
  template <typename T>
  Eigen::Vector3d CalcWeightedPos(const std::vector<double>& weights,
                                  const T& atoms) const;
 
  void PlaceMapAtomonMD(const std::vector<mapAtom*>& fragment_mapatoms,
                        const std::vector<const Atom*>& fragment_mdatoms) const;
 
  void MapMapAtomonMD(const FragInfo& frag,
                      const std::vector<mapAtom*>& fragment_mapatoms,
                      const std::vector<const Atom*>& fragment_mdatoms) const;
 
  // Default, no-op implementation for AtomContainer types whose own
  // mapAtom type has no hasExternalBond()/setExternalBondDirection()
  // interface at all (StaticSegment/PolarSegment's own atom types --
  // this feature is currently QM-specific only). Given an already-
  // computed rotation (rot, no translation -- a direction is
  // translation-invariant by construction, matching the same
  // "rotate the position-minus-reference-point difference only"
  // principle Atom::Rotate/QMAtom::Rotate's own, existing
  // implementation already uses), transfers the MD atom's own,
  // already-recorded, raw external-bond direction onto the mapped
  // atom, rotated into the mapped atom's own reference frame -- and,
  // alongside it, the SEGMENT id (Atom::
  // getExternalBondPartnerSegmentId()) that bond crosses into, copied
  // straight across, unchanged (a segment id needs no rigid-body
  // transform at all, it is not a geometric quantity). Explicitly
  // specialized for SegmentMapper<QMMolecule> below, where QMAtom
  // actually does have this interface -- matching the same,
  // already-established pattern this class already uses for
  // FillMap()/getRank() above.
  void TransferExternalBondDirection(mapAtom*, const Atom*,
                                     const Eigen::Matrix3d&) const {}
 
  // Same default, no-op/specialization pattern as
  // TransferExternalBondDirection above, for the same reason
  // (StaticSegment/PolarSegment's own atom types have no
  // getBondedPartnerIds()/AddBondedPartner() interface at all).
  // Given a "MD-level atom ID -> mapped atom" lookup covering the
  // WHOLE segment (not just this one fragment -- a bonded partner can
  // genuinely be in a different fragment of the same segment), and
  // the corresponding MD-level Atom whose own, raw, MD-level partner
  // IDs need translating, looks each one up in that lookup and, if
  // found (i.e. the partner is itself part of this same segment --
  // if not found, the partner belongs to a genuinely different
  // segment, and is silently skipped, since that partner has no
  // corresponding mapped atom within this call's own Result at all,
  // matching this class's own, existing convention of never crossing
  // segment boundaries within a single map() call), adds its own
  // mapped ID onto map_atom directly. Explicitly specialized for
  // SegmentMapper<QMMolecule> below.
  void TransferBondedPartners(mapAtom*, const Atom*,
                              const std::map<Index, mapAtom*>&) const {}
 
  Logger& log_;
  std::pair<Index, Index> CalcAtomIdRange(const Segment& seg) const;
  std::pair<Index, Index> CalcAtomIdRange(const std::vector<Index>& seg) const;
 
  atom_id StringToMapIndex(const std::string& map_string) const;
 
  atom_id StringToMDIndex(const std::string& md_string) const;
 
  Index getRank(const mapAtom& atom) const { return atom.getRank(); }
 
  std::vector<double> getWeights(const tools::Property& frag) const;
 
  std::string getFrame(const tools::Property& frag) const {
    if (frag.exists(mapatom_xml_.at("frame"))) {
      return frag.get(mapatom_xml_.at("frame")).template as<std::string>();
    }
    return frag.get("localframe").template as<std::string>();
  }
 
  void FillMap() {
    mapatom_xml_["tag"] = "MP";
    mapatom_xml_["name"] = "MPole";
    mapatom_xml_["atoms"] = "mpoles";
    mapatom_xml_["coords"] = "multipoles";
    mapatom_xml_["weights"] = "mp_weights";
    mapatom_xml_["frame"] = "mp_localframe";
  }
};
 
template <>
inline void SegmentMapper<QMMolecule>::FillMap() {
  mapatom_xml_["tag"] = "QM";
  mapatom_xml_["name"] = "QMAtom";
  mapatom_xml_["atoms"] = "qmatoms";
  mapatom_xml_["coords"] = "qmcoords";
  mapatom_xml_["weights"] = "qm_weights";
  mapatom_xml_["frame"] = "qm_localframe";
}
 
template <>
inline Index SegmentMapper<QMMolecule>::getRank(const QMAtom&) const {
  return 0;
}
 
template <>
inline void SegmentMapper<QMMolecule>::TransferExternalBondDirection(
    QMAtom* map_atom, const Atom* md_atom, const Eigen::Matrix3d& rot) const {
  if (md_atom->hasExternalBond()) {
    map_atom->setExternalBondDirection(rot *
                                       md_atom->getExternalBondDirection());
    map_atom->setExternalBondPartnerSegmentId(
        md_atom->getExternalBondPartnerSegmentId());
  }
}
 
template <>
inline void SegmentMapper<QMMolecule>::TransferBondedPartners(
    QMAtom* map_atom, const Atom* md_atom,
    const std::map<Index, QMAtom*>& md_id_to_map_atom) const {
  for (Index i = 0; i < Atom::kMaxBondedPartners; i++) {
    Index partner_md_id = md_atom->getBondedPartnerIds()[i];
    if (partner_md_id == -1) {
      continue;
    }
    auto it = md_id_to_map_atom.find(partner_md_id);
    if (it != md_id_to_map_atom.end()) {
      map_atom->AddBondedPartner(it->second->getId());
    }
  }
}
 
using QMMapper = SegmentMapper<QMMolecule>;
using StaticMapper = SegmentMapper<StaticSegment>;
using PolarMapper = SegmentMapper<PolarSegment>;
}  // namespace xtp
}  // namespace votca
 
#endif  // VOTCA_XTP_SEGMENTMAPPER_H