diff --git a/include/votca/csg/basebead.h b/include/votca/csg/basebead.h index aaa06907aa..e58a67a0f5 100644 --- a/include/votca/csg/basebead.h +++ b/include/votca/csg/basebead.h @@ -14,11 +14,10 @@ * limitations under the License. * */ -#pragma once #ifndef VOTCA_CSG_BASEBEAD_H #define VOTCA_CSG_BASEBEAD_H +#pragma once -#include "topologyitem.h" #include #include #include @@ -74,9 +73,6 @@ class BaseBead { */ Index getMoleculeId() const noexcept { return molecule_id_; } - /// Gets the topology pointer the bead is attached too - Topology *getParent() const { return topology_item_.getParent(); } - /** * get the bead type * \return const string @@ -141,8 +137,6 @@ class BaseBead { protected: BaseBead(){}; - TopologyItem topology_item_ = nullptr; - std::string type_ = tools::topology_constants::unassigned_bead_type; Index id_ = tools::topology_constants::unassigned_residue_id; Index molecule_id_ = tools::topology_constants::unassigned_molecule_id; diff --git a/include/votca/csg/bead.h b/include/votca/csg/bead.h index ac2553b021..6b7ea2853d 100644 --- a/include/votca/csg/bead.h +++ b/include/votca/csg/bead.h @@ -15,8 +15,9 @@ * */ -#ifndef _VOTCA_CSG_BEAD_H -#define _VOTCA_CSG_BEAD_H +#ifndef VOTCA_CSG_BEAD_H +#define VOTCA_CSG_BEAD_H +#pragma once #include #include @@ -289,10 +290,9 @@ class Bead : public BaseBead { bool bead_force_set_; /// constructor - Bead(Topology *owner, Index id, std::string type, Symmetry symmetry, - std::string name, Index resnr, double m, double q) + Bead(Index id, std::string type, Symmetry symmetry, std::string name, + Index resnr, double m, double q) : symmetry_(symmetry), charge_(q), residue_number_(resnr) { - topology_item_._parent = owner; setId(id); setType(type); setName(name); @@ -372,4 +372,4 @@ inline void Bead::HasW(bool b) { bW_ = b; } } // namespace csg } // namespace votca -#endif // _VOTCA_CSG_BEAD_H +#endif // VOTCA_CSG_BEAD_H diff --git a/include/votca/csg/beadtype.h b/include/votca/csg/beadtype.h deleted file mode 100644 index eb08732ee5..0000000000 --- a/include/votca/csg/beadtype.h +++ /dev/null @@ -1,51 +0,0 @@ -/* - * 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. - * - */ - -#ifndef _VOTCA_CSG_BEADTYPE_H -#define _VOTCA_CSG_BEADTYPE_H - -#include "topologyitem.h" -#include - -namespace votca { -namespace csg { - -/** - \brief Bead Type informaton - - Each bead has a type. While the bead name should be unique, - several beads can share the same type. - */ -class BeadType : public TopologyItem { - public: - const Index &getId() const { return _id; } - const std::string &getName() const { return _name; } - void setName(const std::string &name) { _name = name; } - - private: - Index _id; - std::string _name; - - BeadType(Topology *parent, Index id, const std::string &name) - : TopologyItem(parent), _id(id), _name(name) {} - friend class Topology; -}; - -} // namespace csg -} // namespace votca - -#endif /* _VOTCA_CSG_BEADTYPE_H */ diff --git a/include/votca/csg/boundarycondition.h b/include/votca/csg/boundarycondition.h index d2c78a439b..3a06bb91e5 100644 --- a/include/votca/csg/boundarycondition.h +++ b/include/votca/csg/boundarycondition.h @@ -14,9 +14,9 @@ * limitations under the License. * */ -#pragma once #ifndef VOTCA_CSG_BOUNDARYCONDITION_H #define VOTCA_CSG_BOUNDARYCONDITION_H +#pragma once #include #include @@ -62,6 +62,14 @@ class BoundaryCondition { */ const Eigen::Matrix3d &getBox() const noexcept { return _box; }; + /** + * @brief Self explanatory gets the shortest dimension of the boundary + * conditions + * + * @return + */ + double getShortestBoxDimension() const; + /** * get the volume of the box * \return box volume as double diff --git a/include/votca/csg/interaction.h b/include/votca/csg/interaction.h index 26a52dae7d..4a99a2349c 100644 --- a/include/votca/csg/interaction.h +++ b/include/votca/csg/interaction.h @@ -15,8 +15,9 @@ * */ -#ifndef _VOTCA_CSG_INTERACTION_H -#define _VOTCA_CSG_INTERACTION_H +#pragma once +#ifndef VOTCA_CSG_INTERACTION_H +#define VOTCA_CSG_INTERACTION_H #include "bead.h" #include "topology.h" @@ -79,8 +80,8 @@ class Interaction { } virtual Eigen::Vector3d Grad(const Topology &top, Index bead) = 0; - Index BeadCount() { return _beads.size(); } - Index getBeadId(Index bead) { + Index BeadCount() const { return _beads.size(); } + Index getBeadId(Index bead) const { assert(bead > -1 && boost::lexical_cast(bead) < _beads.size()); return _beads[bead]; } @@ -321,4 +322,4 @@ inline Eigen::Vector3d IDihedral::Grad(const Topology &top, Index bead) { } // namespace csg } // namespace votca -#endif // _VOTCA_CSG_INTERACTION_H +#endif // VOTCA_CSG_INTERACTION_H diff --git a/include/votca/csg/map.h b/include/votca/csg/map.h index 8405f8322d..e52539785b 100644 --- a/include/votca/csg/map.h +++ b/include/votca/csg/map.h @@ -15,9 +15,11 @@ * */ -#ifndef _VOTCA_CSG_MAP_H -#define _VOTCA_CSG_MAP_H +#ifndef VOTCA_CSG_MAP_H +#define VOTCA_CSG_MAP_H +#pragma once +#include "boundarycondition.h" #include "molecule.h" #include #include @@ -37,7 +39,7 @@ class Map { void AddBeadMap(BeadMap *bmap) { _maps.push_back(bmap); } - void Apply(); + void Apply(const BoundaryCondition &bc); protected: Molecule _in, _out; @@ -50,7 +52,7 @@ class Map { class BeadMap { public: virtual ~BeadMap() = default; - virtual void Apply() = 0; + virtual void Apply(const BoundaryCondition &) = 0; virtual void Initialize(Molecule *in, Bead *out, tools::Property *opts_map, tools::Property *opts_bead); @@ -76,7 +78,7 @@ inline void BeadMap::Initialize(Molecule *in, Bead *out, class Map_Sphere : public BeadMap { public: Map_Sphere() = default; - void Apply() override; + void Apply(const BoundaryCondition &) override; void Initialize(Molecule *in, Bead *out, tools::Property *opts_bead, tools::Property *opts_map) override; @@ -106,7 +108,7 @@ inline void Map_Sphere::AddElem(Bead *in, double weight, double force_weight) { class Map_Ellipsoid : public Map_Sphere { public: Map_Ellipsoid() = default; - void Apply() override; + void Apply(const BoundaryCondition &) override; protected: }; @@ -114,4 +116,4 @@ class Map_Ellipsoid : public Map_Sphere { } // namespace csg } // namespace votca -#endif /* _VOTCA_CSG_MAP_H */ +#endif // VOTCA_CSG_MAP_H diff --git a/include/votca/csg/molecule.h b/include/votca/csg/molecule.h index 5c34590cda..1cad3994ce 100644 --- a/include/votca/csg/molecule.h +++ b/include/votca/csg/molecule.h @@ -15,11 +15,11 @@ * */ -#ifndef _VOTCA_CSG_MOLECULE_H -#define _VOTCA_CSG_MOLECULE_H +#pragma once +#ifndef VOTCA_CSG_MOLECULE_H +#define VOTCA_CSG_MOLECULE_H #include "bead.h" -#include "topologyitem.h" #include #include #include @@ -39,7 +39,7 @@ class Interaction; \todo sort atoms in molecule */ -class Molecule : public TopologyItem { +class Molecule { public: /// get the molecule ID Index getId() const { return _id; } @@ -98,8 +98,7 @@ class Molecule : public TopologyItem { void *_userdata; /// constructor - Molecule(Topology *parent, Index id, std::string name) - : TopologyItem(parent), _id(id), _name(name) {} + Molecule(Index id, std::string name) : _id(id), _name(name) {} friend class Topology; }; @@ -117,4 +116,4 @@ inline Index Molecule::getBeadIdByName(const std::string &name) { } // namespace csg } // namespace votca -#endif /* _VOTCA_CSG_MOLECULE_H */ +#endif // VOTCA_CSG_MOLECULE_H diff --git a/include/votca/csg/residue.h b/include/votca/csg/residue.h index 84751d693a..ea3ded4f45 100644 --- a/include/votca/csg/residue.h +++ b/include/votca/csg/residue.h @@ -14,11 +14,10 @@ * limitations under the License. * */ +#pragma once +#ifndef VOTCA_CSG_RESIDUE_H +#define VOTCA_CSG_RESIDUE_H -#ifndef _VOTCA_CSG_RESIDUE_H -#define _VOTCA_CSG_RESIDUE_H - -#include "topologyitem.h" #include namespace votca { @@ -32,7 +31,7 @@ namespace csg { based on their residue. */ -class Residue : public TopologyItem { +class Residue { public: /// get the name of the residue const std::string &getName(); @@ -46,8 +45,7 @@ class Residue : public TopologyItem { private: /// constructor - Residue(Topology *parent, Index id, const std::string &name) - : TopologyItem(parent), _id(id), _name(name) {} + Residue(Index id, const std::string &name) : _id(id), _name(name) {} friend class Topology; }; @@ -56,4 +54,4 @@ inline const std::string &Residue::getName() { return _name; } } // namespace csg } // namespace votca -#endif /* _VOTCA_CSG_RESIDUE_H */ +#endif // VOTCA_CSG_RESIDUE_H diff --git a/include/votca/csg/topology.h b/include/votca/csg/topology.h index a4a0d3534b..bbd84efdff 100644 --- a/include/votca/csg/topology.h +++ b/include/votca/csg/topology.h @@ -15,8 +15,9 @@ * */ -#ifndef _VOTCA_CSG_TOPOLOGY_H -#define _VOTCA_CSG_TOPOLOGY_H +#ifndef VOTCA_CSG_TOPOLOGY_H +#define VOTCA_CSG_TOPOLOGY_H +#pragma once #include #include @@ -25,7 +26,6 @@ #include #include "bead.h" -#include "beadtype.h" #include "boundarycondition.h" #include "exclusionlist.h" #include "molecule.h" @@ -176,6 +176,9 @@ class Topology { * @return bonded interaction container */ InteractionContainer &BondedInteractions() { return _interactions; } + const InteractionContainer &BondedInteractions() const { + return _interactions; + } void AddBondedInteraction(Interaction *ic); std::list InteractionsInGroup(const std::string &group); @@ -289,6 +292,13 @@ class Topology { */ const Eigen::Matrix3d &getBox() const { return _bc->getBox(); }; + /** + * @brief Return the boundary condition object + */ + const BoundaryCondition &getBoundary() const { + assert(_bc != nullptr && "Cannot return boundary condition is null"); + return *_bc; + }; /** * set the time of current frame * \param t simulation time in ns @@ -427,25 +437,25 @@ inline Bead *Topology::CreateBead(Bead::Symmetry symmetry, std::string name, std::string type, Index resnr, double m, double q) { - Bead *b = new Bead(this, _beads.size(), type, symmetry, name, resnr, m, q); + Bead *b = new Bead(_beads.size(), type, symmetry, name, resnr, m, q); _beads.push_back(b); return b; } inline Molecule *Topology::CreateMolecule(std::string name) { - Molecule *mol = new Molecule(this, _molecules.size(), name); + Molecule *mol = new Molecule(_molecules.size(), name); _molecules.push_back(mol); return mol; } inline Residue *Topology::CreateResidue(std::string name, Index id) { - Residue *res = new Residue(this, id, name); + Residue *res = new Residue(id, name); _residues.push_back(res); return res; } inline Residue *Topology::CreateResidue(std::string name) { - Residue *res = new Residue(this, _residues.size(), name); + Residue *res = new Residue(_residues.size(), name); _residues.push_back(res); return res; } @@ -464,4 +474,4 @@ inline void Topology::InsertExclusion(Bead *bead1, iteratable &l) { #include "interaction.h" -#endif /* _VOTCA_CSG_TOPOLOGY_H */ +#endif // VOTCA_CSG_TOPOLOGY_H diff --git a/include/votca/csg/topologyitem.h b/include/votca/csg/topologyitem.h deleted file mode 100644 index ef7f18eb91..0000000000 --- a/include/votca/csg/topologyitem.h +++ /dev/null @@ -1,44 +0,0 @@ -/* - * 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. - * - */ - -#ifndef _VOTCA_CSG_TOPOLOGYITEM_H -#define _VOTCA_CSG_TOPOLOGYITEM_H - -namespace votca { -namespace csg { - -class Topology; - -class TopologyItem { - public: - virtual ~TopologyItem() = default; - Topology *getParent() const { return _parent; } - - protected: - TopologyItem(Topology *parent) : _parent(parent) {} - - Topology *_parent; - - friend class Topology; - friend class BaseBead; - friend class Bead; -}; - -} // namespace csg -} // namespace votca - -#endif /* _VOTCA_CSG_TOPOLOGYITEM_H */ diff --git a/src/csg_boltzmann/main.cc b/src/csg_boltzmann/main.cc index 06c41e611b..751482ebb3 100644 --- a/src/csg_boltzmann/main.cc +++ b/src/csg_boltzmann/main.cc @@ -53,7 +53,9 @@ class CsgBoltzmann : public CsgApplication { bool EvaluateTopology(Topology *top, Topology *top_ref) override; protected: - ExclusionList *CreateExclusionList(Molecule &atomistic, Molecule &cg); + ExclusionList *CreateExclusionList(Topology *top_atomistic, + Molecule &atomistic, Topology *top_cg, + Molecule &cg); BondedStatistics _bs; }; void CsgBoltzmann::Initialize() { @@ -86,7 +88,7 @@ bool CsgBoltzmann::EvaluateTopology(Topology *top, Topology *top_ref) { << " in coarse grained representation " << top_ref->MoleculeByIndex(0)->getName() << endl; - ex = CreateExclusionList(*top_ref->MoleculeByIndex(0), + ex = CreateExclusionList(top_ref, *top_ref->MoleculeByIndex(0), top, *top->MoleculeByIndex(0)); std::ofstream fl; fl.open(OptionsMap()["excl"].as()); @@ -101,7 +103,9 @@ bool CsgBoltzmann::EvaluateTopology(Topology *top, Topology *top_ref) { return true; } -ExclusionList *CsgBoltzmann::CreateExclusionList(Molecule &atomistic, +ExclusionList *CsgBoltzmann::CreateExclusionList(Topology *top_atomistic, + Molecule &atomistic, + Topology *top_cg, Molecule &cg) { ExclusionList *ex = new ExclusionList(); // exclude all with all @@ -114,23 +118,21 @@ ExclusionList *CsgBoltzmann::CreateExclusionList(Molecule &atomistic, } // remove exclusions from inside a mapped bead - Topology *at_top = atomistic.getParent(); for (votca::Index i = 0; i < cg.BeadCount(); ++i) { const vector &parent_beads = cg.getBead(i)->ParentBeads(); list excl_list; for (const votca::Index &parent_bead_id : parent_beads) { - excl_list.push_back(at_top->getBead(parent_bead_id)); + excl_list.push_back(top_atomistic->getBead(parent_bead_id)); } ex->Remove(excl_list); } // remove exclusion which come from atomistic topology and hence bonds and // angles - Topology *cg_top = cg.getParent(); for (votca::Index i = 0; i < cg.BeadCount() - 1; ++i) { for (votca::Index j = i + 1; j < cg.BeadCount(); ++j) { - if (cg_top->getExclusions().IsExcluded(cg.getBead(i), cg.getBead(j))) { + if (top_cg->getExclusions().IsExcluded(cg.getBead(i), cg.getBead(j))) { const vector &parent_beads_w = cg.getBead(i)->ParentBeads(); const vector &parent_beads_v = @@ -138,8 +140,8 @@ ExclusionList *CsgBoltzmann::CreateExclusionList(Molecule &atomistic, for (const votca::Index parent_bead_id_w : parent_beads_w) { for (const votca::Index parent_bead_id_v : parent_beads_v) { - ex->RemoveExclusion(at_top->getBead(parent_bead_id_w), - at_top->getBead(parent_bead_id_v)); + ex->RemoveExclusion(top_atomistic->getBead(parent_bead_id_w), + top_atomistic->getBead(parent_bead_id_v)); } } } diff --git a/src/libcsg/boundarycondition.cc b/src/libcsg/boundarycondition.cc index 79cf529bd0..32fd97628e 100644 --- a/src/libcsg/boundarycondition.cc +++ b/src/libcsg/boundarycondition.cc @@ -14,9 +14,10 @@ * limitations under the License. * */ +#include +#include #include "../../include/votca/csg/boundarycondition.h" -#include namespace votca { namespace csg { @@ -25,5 +26,29 @@ double BoundaryCondition::BoxVolume() const noexcept { return std::abs(_box.determinant()); } +double BoundaryCondition::getShortestBoxDimension() const { + assert(getBoxType() != eBoxtype::typeOpen && + "Cannot get the shortest dimension of the box because it is open"); + + Eigen::Vector3d box_a = _box.col(0); + Eigen::Vector3d box_b = _box.col(1); + Eigen::Vector3d box_c = _box.col(2); + + // create plane normals + Eigen::Vector3d norm_a = box_b.cross(box_c); + Eigen::Vector3d norm_b = box_c.cross(box_a); + Eigen::Vector3d norm_c = box_a.cross(box_b); + + norm_a.normalize(); + norm_b.normalize(); + norm_c.normalize(); + + double la = box_a.dot(norm_a); + double lb = box_b.dot(norm_b); + double lc = box_c.dot(norm_c); + + return std::min(la, std::min(lb, lc)); +} + } // namespace csg } // namespace votca diff --git a/src/libcsg/csgapplication.cc b/src/libcsg/csgapplication.cc index 79fdc2585c..d279da33a3 100644 --- a/src/libcsg/csgapplication.cc +++ b/src/libcsg/csgapplication.cc @@ -225,6 +225,9 @@ void CsgApplication::Run(void) { // read in the topology for master ////////////////////////////////////////////////// reader->ReadTopology(_op_vm["top"].as(), master->_top); + // Ensure that the coarse grained topology will have the same boundaries + master->_top_cg.setBox(master->_top.getBox()); + std::cout << "I have " << master->_top.BeadCount() << " beads in " << master->_top.MoleculeCount() << " molecules" << std::endl; master->_top.CheckMoleculeNaming(); @@ -249,9 +252,17 @@ void CsgApplication::Run(void) { std::cout << "I have " << master->_top_cg.BeadCount() << " beads in " << master->_top_cg.MoleculeCount() << " molecules for the coarsegraining" << std::endl; - master->_map->Apply(); - if (!EvaluateTopology(&master->_top_cg, &master->_top)) { - return; + + // If the trajectory reader is off but mapping flag is specified do apply + // the mapping, this switch is necessary in cases where xml files are + // specified, which do not contain positional information. In such cases + // it is not possible to apply the positional mapping, a trajectory file + // must be read in. + if (DoTrajectory() == false) { + master->_map->Apply(); + if (!EvaluateTopology(&master->_top_cg, &master->_top)) { + return; + } } } else if (!EvaluateTopology(&master->_top)) { return; diff --git a/src/libcsg/map.cc b/src/libcsg/map.cc index 2c6605501f..96edc3925a 100644 --- a/src/libcsg/map.cc +++ b/src/libcsg/map.cc @@ -45,9 +45,9 @@ Map::~Map() { _maps.clear(); } -void Map::Apply() { +void Map::Apply(const BoundaryCondition &bc) { for (auto &_map : _maps) { - _map->Apply(); + _map->Apply(bc); } } @@ -129,15 +129,15 @@ void Map_Sphere::Initialize(Molecule *in, Bead *out, Property *opts_bead, } } -void Map_Sphere::Apply() { +void Map_Sphere::Apply(const BoundaryCondition &bc) { + + assert(_matrix.size() > 0 && "Cannot map to sphere there are no beads"); bool bPos, bVel, bF; bPos = bVel = bF = false; _out->ClearParentBeads(); // the following is needed for pbc treatment - Topology *top = _out->getParent(); - double max_dist = 0.5 * top->ShortestBoxSize(); Eigen::Vector3d r0 = Eigen::Vector3d::Zero(); string name0; Index id0 = 0; @@ -154,25 +154,45 @@ void Map_Sphere::Apply() { Eigen::Vector3d f = Eigen::Vector3d::Zero(); Eigen::Vector3d vel = Eigen::Vector3d::Zero(); + Bead *bead_max_dist = _matrix.at(0)._in; + double max_bead_dist = + bc.BCShortestConnection(r0, bead_max_dist->getPos()).norm(); + for (auto &iter : _matrix) { Bead *bead = iter._in; _out->AddParentBead(bead->getId()); M += bead->getMass(); if (bead->HasPos()) { - Eigen::Vector3d r = top->BCShortestConnection(r0, bead->getPos()); - if (r.norm() > max_dist) { - cout << r0 << " " << bead->getPos() << endl; - throw std::runtime_error( - "coarse-grained bead is bigger than half the box \n (atoms " + - name0 + " (id " + boost::lexical_cast(id0 + 1) + ")" + - ", " + bead->getName() + " (id " + - boost::lexical_cast(bead->getId() + 1) + ")" + - +" , molecule " + - boost::lexical_cast(bead->getMoleculeId() + 1) + ")"); + Eigen::Vector3d r = bc.BCShortestConnection(r0, bead->getPos()); + if (r.norm() > max_bead_dist) { + max_bead_dist = r.norm(); + bead_max_dist = bead; } cg += iter._weight * (r + r0); bPos = true; } + } + + /// Safety check, if box is not open check if the bead is larger than the + /// boundaries + if (bc.getBoxType() != BoundaryCondition::eBoxtype::typeOpen) { + double max_dist = 0.5 * bc.getShortestBoxDimension(); + if (max_bead_dist > max_dist) { + cout << r0 << " " << bead_max_dist->getPos() << endl; + throw std::runtime_error( + "coarse-grained bead is bigger than half the box \n " + "(atoms " + + name0 + " (id " + boost::lexical_cast(id0 + 1) + ")" + ", " + + bead_max_dist->getName() + " (id " + + boost::lexical_cast(bead_max_dist->getId() + 1) + ")" + + +" , molecule " + + boost::lexical_cast(bead_max_dist->getMoleculeId() + 1) + + ")"); + } + } + + for (auto &iter : _matrix) { + Bead *bead = iter._in; if (bead->HasVel()) { vel += iter._weight * bead->getVel(); bVel = true; @@ -195,18 +215,22 @@ void Map_Sphere::Apply() { } /// \todo implement this function -void Map_Ellipsoid::Apply() { +void Map_Ellipsoid::Apply(const BoundaryCondition &bc) { + + assert(_matrix.size() > 0 && "Cannot map to ellipsoid there are no beads"); bool bPos, bVel, bF; bPos = bVel = bF = false; // the following is needed for pbc treatment - Topology *top = _out->getParent(); - double max_dist = 0.5 * top->ShortestBoxSize(); Eigen::Vector3d r0 = Eigen::Vector3d::Zero(); + string name0; + Index id0 = 0; if (_matrix.size() > 0) { if (_matrix.front()._in->HasPos()) { r0 = _matrix.front()._in->getPos(); + name0 = _matrix.front()._in->getName(); + id0 = _matrix.front()._in->getId(); } } Eigen::Vector3d cg = Eigen::Vector3d::Zero(); @@ -217,18 +241,45 @@ void Map_Ellipsoid::Apply() { Index n; n = 0; _out->ClearParentBeads(); + + Bead *bead_max_dist = _matrix.at(0)._in; + double max_bead_dist = + bc.BCShortestConnection(r0, bead_max_dist->getPos()).norm(); + for (auto &iter : _matrix) { Bead *bead = iter._in; _out->AddParentBead(bead->getId()); if (bead->HasPos()) { - Eigen::Vector3d r = top->BCShortestConnection(r0, bead->getPos()); - if (r.norm() > max_dist) { - throw std::runtime_error( - "coarse-grained bead is bigger than half the box"); + Eigen::Vector3d r = bc.BCShortestConnection(r0, bead->getPos()); + if (r.norm() > max_bead_dist) { + max_bead_dist = r.norm(); + bead_max_dist = bead; } cg += iter._weight * (r + r0); bPos = true; } + } + + /// Safety check, if box is not open check if the bead is larger than the + /// boundaries + if (bc.getBoxType() != BoundaryCondition::eBoxtype::typeOpen) { + double max_dist = 0.5 * bc.getShortestBoxDimension(); + if (max_bead_dist > max_dist) { + cout << r0 << " " << bead_max_dist->getPos() << endl; + throw std::runtime_error( + "coarse-grained bead is bigger than half the box \n " + "(atoms " + + name0 + " (id " + boost::lexical_cast(id0 + 1) + ")" + ", " + + bead_max_dist->getName() + " (id " + + boost::lexical_cast(bead_max_dist->getId() + 1) + ")" + + +" , molecule " + + boost::lexical_cast(bead_max_dist->getMoleculeId() + 1) + + ")"); + } + } + + for (auto &iter : _matrix) { + Bead *bead = iter._in; if (bead->HasVel() == true) { vel += iter._weight * bead->getVel(); bVel = true; diff --git a/src/libcsg/topology.cc b/src/libcsg/topology.cc index 96bd6f7ae3..b6e0abe4ed 100644 --- a/src/libcsg/topology.cc +++ b/src/libcsg/topology.cc @@ -336,24 +336,7 @@ BoundaryCondition::eBoxtype Topology::autoDetectBoxType( } double Topology::ShortestBoxSize() const { - Eigen::Vector3d box_a = getBox().col(0); - Eigen::Vector3d box_b = getBox().col(1); - Eigen::Vector3d box_c = getBox().col(2); - - // create plane normals - Eigen::Vector3d norm_a = box_b.cross(box_c); - Eigen::Vector3d norm_b = box_c.cross(box_a); - Eigen::Vector3d norm_c = box_a.cross(box_b); - - norm_a.normalize(); - norm_b.normalize(); - norm_c.normalize(); - - double la = box_a.dot(norm_a); - double lb = box_b.dot(norm_b); - double lc = box_c.dot(norm_c); - - return std::min(la, std::min(lb, lc)); + return _bc->getShortestBoxDimension(); } } // namespace csg diff --git a/src/libcsg/topologymap.cc b/src/libcsg/topologymap.cc index 1c55f60ad0..a707e08eca 100644 --- a/src/libcsg/topologymap.cc +++ b/src/libcsg/topologymap.cc @@ -16,6 +16,7 @@ */ #include "../../include/votca/csg/topologymap.h" +#include "../../include/votca/csg/boundarycondition.h" namespace votca { namespace csg { @@ -34,7 +35,7 @@ void TopologyMap::Apply() { _out->setBox(_in->getBox()); for (auto& _map : _maps) { - _map->Apply(); + _map->Apply(_out->getBoundary()); } } diff --git a/src/tests/test_basebead.cc b/src/tests/test_basebead.cc index c617649d58..cf0617c297 100644 --- a/src/tests/test_basebead.cc +++ b/src/tests/test_basebead.cc @@ -19,7 +19,6 @@ #define BOOST_TEST_MODULE basebead_test #include "../../include/votca/csg/basebead.h" -#include "../../include/votca/csg/beadtype.h" #include "../../include/votca/csg/molecule.h" #include "../../include/votca/csg/topology.h" #include diff --git a/src/tests/test_bead.cc b/src/tests/test_bead.cc index a0fb9e79d5..42a0bdfdc4 100644 --- a/src/tests/test_bead.cc +++ b/src/tests/test_bead.cc @@ -19,7 +19,6 @@ #define BOOST_TEST_MODULE bead_test #include "../../include/votca/csg/bead.h" -#include "../../include/votca/csg/beadtype.h" #include "../../include/votca/csg/molecule.h" #include "../../include/votca/csg/topology.h" #include diff --git a/src/tests/test_beadtriple.cc b/src/tests/test_beadtriple.cc index 0303dbccef..f24a3c79f9 100644 --- a/src/tests/test_beadtriple.cc +++ b/src/tests/test_beadtriple.cc @@ -22,7 +22,6 @@ #include "../../include/votca/csg/bead.h" #include "../../include/votca/csg/beadtriple.h" -#include "../../include/votca/csg/beadtype.h" #include "../../include/votca/csg/topology.h" #include #include diff --git a/src/tests/test_interaction.cc b/src/tests/test_interaction.cc index 735a0c83cf..9ca62902be 100644 --- a/src/tests/test_interaction.cc +++ b/src/tests/test_interaction.cc @@ -21,7 +21,6 @@ #include #include "../../include/votca/csg/bead.h" -#include "../../include/votca/csg/beadtype.h" #include "../../include/votca/csg/interaction.h" #include "../../include/votca/csg/molecule.h" #include "../../include/votca/csg/topology.h" diff --git a/src/tests/test_nblist_3body.cc b/src/tests/test_nblist_3body.cc index 33bac14518..c4f07c085f 100644 --- a/src/tests/test_nblist_3body.cc +++ b/src/tests/test_nblist_3body.cc @@ -22,7 +22,6 @@ #include "../../include/votca/csg/bead.h" #include "../../include/votca/csg/beadlist.h" -#include "../../include/votca/csg/beadtype.h" #include "../../include/votca/csg/nblist_3body.h" #include "../../include/votca/csg/topology.h" #include diff --git a/src/tests/test_nblistgrid_3body.cc b/src/tests/test_nblistgrid_3body.cc index 86606b59c3..46a8d0fa6f 100644 --- a/src/tests/test_nblistgrid_3body.cc +++ b/src/tests/test_nblistgrid_3body.cc @@ -22,7 +22,6 @@ #include "../../include/votca/csg/bead.h" #include "../../include/votca/csg/beadlist.h" -#include "../../include/votca/csg/beadtype.h" #include "../../include/votca/csg/nblistgrid_3body.h" #include "../../include/votca/csg/topology.h" #include diff --git a/src/tests/test_triplelist.cc b/src/tests/test_triplelist.cc index 465119531b..5ff9039b65 100644 --- a/src/tests/test_triplelist.cc +++ b/src/tests/test_triplelist.cc @@ -20,7 +20,6 @@ #define BOOST_TEST_MODULE triplelist_test #include "../../include/votca/csg/bead.h" #include "../../include/votca/csg/beadtriple.h" -#include "../../include/votca/csg/beadtype.h" #include "../../include/votca/csg/topology.h" #include "../../include/votca/csg/triplelist.h" #include diff --git a/src/tools/csg_gmxtopol.cc b/src/tools/csg_gmxtopol.cc index 5b7f693a25..9e9e387552 100644 --- a/src/tools/csg_gmxtopol.cc +++ b/src/tools/csg_gmxtopol.cc @@ -44,8 +44,8 @@ class GmxTopolApp : public CsgApplication { protected: void WriteAtoms(ostream &out, Molecule &cg); - void WriteInteractions(ostream &out, Molecule &cg); - void WriteMolecule(ostream &out, Molecule &cg); + void WriteInteractions(ostream &out, const Topology &top, Molecule &cg); + void WriteMolecule(ostream &out, const Topology &top, Molecule &cg); }; void GmxTopolApp::Initialize(void) { @@ -62,7 +62,7 @@ bool GmxTopolApp::EvaluateTopology(Topology *top, Topology *) { } ofstream fl; fl.open((OptionsMap()["out"].as() + ".top")); - WriteMolecule(fl, *(top->MoleculeByIndex(0))); + WriteMolecule(fl, *top, *(top->MoleculeByIndex(0))); fl.close(); return true; } @@ -78,10 +78,11 @@ void GmxTopolApp::WriteAtoms(ostream &out, Molecule &cg) { out << endl; } -void GmxTopolApp::WriteInteractions(ostream &out, Molecule &cg) { +void GmxTopolApp::WriteInteractions(ostream &out, const Topology &top, + Molecule &cg) { votca::Index nb = -1; - for (Interaction *ic : cg.getParent()->BondedInteractions()) { + for (const Interaction *ic : top.BondedInteractions()) { if (ic->getMolecule() != cg.getId()) { continue; } @@ -114,12 +115,13 @@ void GmxTopolApp::WriteInteractions(ostream &out, Molecule &cg) { } } -void GmxTopolApp::WriteMolecule(ostream &out, Molecule &cg) { +void GmxTopolApp::WriteMolecule(ostream &out, const Topology &top, + Molecule &cg) { out << "[ moleculetype ]\n"; out << cg.getName() << " 3\n\n"; WriteAtoms(out, cg); - WriteInteractions(out, cg); + WriteInteractions(out, top, cg); } int main(int argc, char **argv) {