topology.h

Go to the documentation of this file.

/*
 * Copyright 2009-2021 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.
 *
 */
 
#ifndef VOTCA_CSG_TOPOLOGY_H
#define VOTCA_CSG_TOPOLOGY_H
#pragma once
 
// Standard includes
#include <cassert>
#include <map>
#include <memory>
#include <unordered_map>
#include <vector>
 
// Third party includes
#include <boost/container/deque.hpp>
 
// VOTCA includes
#include <votca/tools/types.h>
 
// Local VOTCA includes
#include "bead.h"
#include "boundarycondition.h"
#include "exclusionlist.h"
#include "molecule.h"
#include "openbox.h"
#include "orthorhombicbox.h"
#include "residue.h"
#include "triclinicbox.h"
 
namespace votca {
namespace csg {
 
class Interaction;
 
/* Boost deque has been chosen and contents have been replaced with objects
 * as opposed to heap allocated types:
 * 1. To get rid of indirection
 * 2. Clarify ownership
 * 3. To ensure pointers are not invalidated if the container changes size
 * that is not a guarantee with a vector
 * 4. To provide better contiguous memory access, not possible with std::deque
 * or list
 */
typedef boost::container::deque_options<
    boost::container::block_size<sizeof(Residue) * 4>>::type block_residue_x4_t;
typedef boost::container::deque_options<
    boost::container::block_size<sizeof(Molecule) * 4>>::type
    block_molecule_4x_t;
typedef boost::container::deque_options<
    boost::container::block_size<sizeof(Bead) * 4>>::type block_bead_x4_t;
 
using MoleculeContainer =
    boost::container::deque<Molecule, void, block_molecule_4x_t>;
using BeadContainer = boost::container::deque<Bead, void, block_bead_x4_t>;
using ResidueContainer =
    boost::container::deque<Residue, void, block_residue_x4_t>;
using InteractionContainer = std::vector<Interaction *>;

class Topology {
 public:
  Topology() { bc_ = std::make_unique<OpenBox>(); }
  ~Topology();

  void Cleanup();

  Bead *CreateBead(Bead::Symmetry symmetry, std::string name, std::string type,
                   Index resnr, double m, double q);

  Molecule *CreateMolecule(std::string name);

  void CheckMoleculeNaming(void);

  Residue &CreateResidue(std::string name);
  Residue &CreateResidue(std::string name, Index id);

  void CreateMoleculesByRange(std::string name, Index first, Index nbeads,
                              Index nmolecules);

  Index MoleculeCount() const { return molecules_.size(); }

  Index BeadCount() const { return beads_.size(); }

  Index ResidueCount() const { return residues_.size(); }

  Molecule *MoleculeByIndex(Index index);

  BeadContainer &Beads() { return beads_; }

  ResidueContainer &Residues() { return residues_; }
  const ResidueContainer &Residues() const { return residues_; }

  MoleculeContainer &Molecules() { return molecules_; }
  const MoleculeContainer &Molecules() const { return molecules_; }

  InteractionContainer &BondedInteractions() { return interactions_; }
  const InteractionContainer &BondedInteractions() const {
    return interactions_;
  }
 
  void AddBondedInteraction(Interaction *ic);
  std::vector<Interaction *> InteractionsInGroup(const std::string &group);

  bool BeadTypeExist(std::string type) const;

  void RegisterBeadType(std::string type);

  Index getBeadTypeId(std::string type) const;

  Bead *getBead(const Index i) { return &beads_[i]; }
  const Bead *getBead(const Index i) const { return &beads_[i]; }
  Residue &getResidue(const Index i) { return residues_[i]; }
  const Residue &getResidue(const Index i) const { return residues_[i]; }
  Molecule *getMolecule(const Index i) { return &molecules_[i]; }
  const Molecule *getMolecule(const Index i) const { return &molecules_[i]; }

  void ClearMoleculeList() { molecules_.clear(); }

  void CopyTopologyData(Topology *top);

  void RenameMolecules(std::string range, std::string name);

  void RenameBeadType(std::string name, std::string newname);

  void SetBeadTypeMass(std::string name, double value);

  void setBox(const Eigen::Matrix3d &box, BoundaryCondition::eBoxtype boxtype =
                                              BoundaryCondition::typeAuto) {
    // determine box type automatically in case boxtype==typeAuto
    if (boxtype == BoundaryCondition::typeAuto) {
      boxtype = autoDetectBoxType(box);
    }
 
    switch (boxtype) {
      case BoundaryCondition::typeTriclinic:
        bc_ = std::make_unique<TriclinicBox>();
        break;
      case BoundaryCondition::typeOrthorhombic:
        bc_ = std::make_unique<OrthorhombicBox>();
        break;
      default:
        bc_ = std::make_unique<OpenBox>();
        break;
    }
 
    bc_->setBox(box);
  };

  const Eigen::Matrix3d &getBox() const { return bc_->getBox(); };

  const BoundaryCondition &getBoundary() const {
    assert(bc_ && "Cannot return boundary condition is null");
    return *bc_;
  };

  void setTime(double t) { time_ = t; };

  double getTime() const { return time_; };

  void setStep(Index s) { step_ = s; };

  Index getStep() const { return step_; };

  void setParticleGroup(std::string particle_group) {
    particle_group_ = particle_group;
  };

  std::string getParticleGroup() const { return particle_group_; };

  Eigen::Vector3d getDist(Index bead1, Index bead2) const;

  Eigen::Vector3d BCShortestConnection(const Eigen::Vector3d &r_i,
                                       const Eigen::Vector3d &r_j) const;

  double ShortestBoxSize() const;

  double BoxVolume() const;

  void RebuildExclusions();

  ExclusionList &getExclusions() { return exclusions_; }
  const ExclusionList &getExclusions() const { return exclusions_; }
 
  BoundaryCondition::eBoxtype getBoxType() const { return bc_->getBoxType(); }
 
  template <typename iteratable>
  void InsertExclusion(Bead *bead1, iteratable &l);
 
  bool HasVel() { return has_vel_; }
  void SetHasVel(const bool v) { has_vel_ = v; }
 
  bool HasForce() { return has_force_; }
  void SetHasForce(const bool v) { has_force_ = v; }
 
 protected:
  std::unique_ptr<BoundaryCondition> bc_;
 
  BoundaryCondition::eBoxtype autoDetectBoxType(
      const Eigen::Matrix3d &box) const;

  std::unordered_map<std::string, Index> beadtypes_;

  BeadContainer beads_;

  MoleculeContainer molecules_;

  ResidueContainer residues_;

  InteractionContainer interactions_;
 
  ExclusionList exclusions_;
 
  std::map<std::string, Index> interaction_groups_;
 
  std::map<std::string, std::vector<Interaction *>> interactions_by_group_;
 
  double time_ = 0.0;
  Index step_ = 0;
  bool has_vel_ = false;
  bool has_force_ = false;

  std::string particle_group_ = "unassigned";
};
 
inline Bead *Topology::CreateBead(Bead::Symmetry symmetry, std::string name,
                                  std::string type, Index resnr, double m,
                                  double q) {
 
  beads_.push_back(Bead(beads_.size(), type, symmetry, name, resnr, m, q));
  return &beads_.back();
}
 
inline Molecule *Topology::CreateMolecule(std::string name) {
  molecules_.push_back(Molecule(molecules_.size(), name));
  return &molecules_.back();
}
 
inline Residue &Topology::CreateResidue(std::string name, Index id) {
  // Note that Residue constructor is intentionally private and only topology
  // class can create it, hence emplace back will not work because the vector
  // class does not have access to the constructor.
  residues_.push_back(Residue(id, name));
  return residues_.back();
}
 
inline Residue &Topology::CreateResidue(std::string name) {
  // Note that Residue constructor is intentionally private and only topology
  // class can create it, hence emplace back will not work because the vector
  // class does not have access to the constructor.
  residues_.push_back(Residue(residues_.size(), name));
  return residues_.back();
}
 
inline Molecule *Topology::MoleculeByIndex(Index index) {
  return &molecules_[index];
}
 
template <typename iteratable>
inline void Topology::InsertExclusion(Bead *bead1, iteratable &l) {
  exclusions_.InsertExclusion(bead1, l);
}
 
}  // namespace csg
}  // namespace votca
 
#include "interaction.h"
 
#endif  // VOTCA_CSG_TOPOLOGY_H