Codebase list votca-xtp / debian/1.6.3-1 include / votca / xtp / atomcontainer.h
debian/1.6.3-1

Tree @debian/1.6.3-1 (Download .tar.gz)

atomcontainer.h @debian/1.6.3-1raw · history · blame

/*
 *            Copyright 2009-2019 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_ATOMCONTAINER_H
#define VOTCA_XTP_ATOMCONTAINER_H
#include <limits>
#include <typeinfo>
#include <votca/tools/elements.h>
#include <votca/xtp/checkpoint.h>
#include <votca/xtp/eigen.h>

/**
 * \brief Basic Container for QMAtoms,PolarSites and Atoms
 *
 *
 *
 */

namespace votca {
namespace xtp {

template <class T>
class AtomContainer {
 public:
  AtomContainer(std::string type, Index id) : _type(type), _id(id){};

  AtomContainer(CheckpointReader& r) { this->ReadFromCpt(r); }
  virtual ~AtomContainer() = default;

  using iterator = typename std::vector<T>::iterator;

  const std::string& getType() const { return _type; }

  void setType(std::string type) { _type = type; }

  Index getId() const { return _id; }

  Index size() const { return _atomlist.size(); }

  void push_back(const T& atom) {
    _atomlist.push_back(atom);
    calcPos();
  }
  void push_back(T&& atom) {
    _atomlist.push_back(atom);
    calcPos();
  }

  void AddContainer(const AtomContainer<T>& container) {
    _type += "_" + container._type;
    _atomlist.insert(_atomlist.end(), container._atomlist.begin(),
                     container._atomlist.end());
    calcPos();
  }

  const T& at(Index index) const { return _atomlist.at(index); }
  T& at(Index index) { return _atomlist.at(index); }

  const T& operator[](Index index) const { return _atomlist[index]; }
  T& operator[](Index index) { return _atomlist[index]; }

  typename std::vector<T>::iterator begin() { return _atomlist.begin(); }
  typename std::vector<T>::iterator end() { return _atomlist.end(); }

  typename std::vector<T>::const_iterator begin() const {
    return _atomlist.begin();
  }
  typename std::vector<T>::const_iterator end() const {
    return _atomlist.end();
  }

  const Eigen::Vector3d& getPos() const { return _pos; }

  // calculates the lowest and highest point in the cube, sorrounding the
  // molecule
  std::pair<Eigen::Vector3d, Eigen::Vector3d> CalcSpatialMinMax() const {
    std::pair<Eigen::Vector3d, Eigen::Vector3d> result;
    Eigen::Vector3d min =
        std::numeric_limits<double>::max() * Eigen::Vector3d::Ones();
    Eigen::Vector3d max =
        std::numeric_limits<double>::min() * Eigen::Vector3d::Ones();
    for (const T& atom : _atomlist) {
      const Eigen::Vector3d& pos = atom.getPos();
      if (pos.x() < min.x()) {
        min.x() = pos.x();
      }
      if (pos.x() > max.x()) {
        max.x() = pos.x();
      }
      if (pos.y() < min.y()) {
        min.y() = pos.y();
      }
      if (pos.y() > max.y()) {
        max.y() = pos.y();
      }
      if (pos.z() < min.z()) {
        min.z() = pos.z();
      }
      if (pos.z() > max.z()) {
        max.z() = pos.z();
      }
    }
    result.first = min;
    result.second = max;
    return result;
  }

  std::vector<std::string> FindUniqueElements() const {
    std::vector<std::string> result;
    for (const T& atom : _atomlist) {
      if (std::find(result.begin(), result.end(), atom.getElement()) ==
          result.end()) {
        result.push_back(atom.getElement());
      }
    }
    return result;
  }

  void Translate(const Eigen::Vector3d& shift) {
    for (T& atom : _atomlist) {
      atom.Translate(shift);
    }
    _pos += shift;
  }

  void Rotate(const Eigen::Matrix3d& R, const Eigen::Vector3d& ref_pos) {
    for (T& atom : _atomlist) {
      atom.Rotate(R, ref_pos);
    }
    calcPos();
  }

  virtual void WriteToCpt(CheckpointWriter& w) const {
    w(_type, "type");
    w(_id, "id");
    w(int(_atomlist.size()), "size");
    T element(0, "H", Eigen::Vector3d::Zero());
    CptTable table =
        w.openTable(element.identify() + "s", element, _atomlist.size());
    std::vector<typename T::data> dataVec(_atomlist.size());
    for (std::size_t i = 0; i < _atomlist.size(); ++i) {
      _atomlist[i].WriteData(dataVec[i]);
    }

    table.write(dataVec);
  }
  virtual void ReadFromCpt(CheckpointReader& r) {
    r(_type, "type");
    r(_id, "id");
    Index size = 0;
    r(size, "size");
    if (size == 0) {
      return;
    }
    T element(0, "H", Eigen::Vector3d::Zero());  // dummy element to get
                                                 // .identify for type
    CptTable table = r.openTable(element.identify() + "s", element);
    _atomlist.clear();
    _atomlist.reserve(table.numRows());
    std::vector<typename T::data> dataVec(table.numRows());
    table.read(dataVec);
    for (std::size_t i = 0; i < table.numRows(); ++i) {
      _atomlist.push_back(T(dataVec[i]));
    }
    calcPos();
  }

  void calcPos() {
    tools::Elements element;
    Eigen::Vector3d pos = Eigen::Vector3d::Zero();
    double totalmass = 0.0;
    for (const T& atom : _atomlist) {
      double mass = element.getMass(atom.getElement());
      totalmass += mass;
      pos += mass * atom.getPos();
    }
    _pos = pos / totalmass;
  }

 protected:
  std::vector<T> _atomlist;
  std::string _type;
  Index _id;

 private:
  Eigen::Vector3d _pos = Eigen::Vector3d::Zero();
};
}  // namespace xtp
}  // namespace votca

#endif  // VOTCA_XTP_ATOMCONTAINER_H