ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
bonded_interaction_data.hpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-2022 The ESPResSo project
3 *
4 * This file is part of ESPResSo.
5 *
6 * ESPResSo is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * ESPResSo is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 */
19
20#pragma once
21
22/** @file
23 * Data structures for bonded interactions.
24 * For more information on how to add new interactions, see @ref bondedIA_new.
25 */
26
27#include "angle_common.hpp"
28#include "angle_cosine.hpp"
29#include "angle_cossquare.hpp"
30#include "angle_harmonic.hpp"
31#include "bonded_coulomb.hpp"
32#include "bonded_coulomb_sr.hpp"
33#include "bonded_tab.hpp"
34#include "dihedral.hpp"
35#include "fene.hpp"
36#include "harmonic.hpp"
42#include "quartic.hpp"
43#include "rigid_bond.hpp"
44#include "thermalized_bond.hpp"
45
46#include "BondList.hpp"
48#include "system/Leaf.hpp"
49
50#include <boost/variant.hpp>
51
52#include <algorithm>
53#include <cassert>
54#include <cmath>
55#include <optional>
56#include <stdexcept>
57#include <unordered_map>
58#include <vector>
59
60/* Special cutoff value for a disabled bond.
61 * Bonds that have this cutoff are not visited during bond evaluation.
62 */
63static constexpr double BONDED_INACTIVE_CUTOFF = -1.;
64
65/** Interaction type for unused bonded interaction slots */
66struct NoneBond {
67 static constexpr int num = 0;
68 double cutoff() const { return BONDED_INACTIVE_CUTOFF; }
69};
70
71/** Interaction type for virtual bonds */
73 static constexpr int num = 1;
74 double cutoff() const { return BONDED_INACTIVE_CUTOFF; }
75};
76
77/** Visitor to get the number of bound partners from the bond parameter
78 * variant.
79 */
80class BondNumPartners : public boost::static_visitor<int> {
81public:
82 template <typename T> int operator()(T const &) const { return T::num; }
83};
84
85/** Variant in which to store the parameters of an individual bonded
86 * interaction
87 */
95
96/**
97 * @brief container for bonded interactions.
98 */
99class BondedInteractionsMap : public System::Leaf<BondedInteractionsMap> {
100 using container_type =
101 std::unordered_map<int, std::shared_ptr<Bonded_IA_Parameters>>;
102
103public:
104 using key_type = typename container_type::key_type;
105 using mapped_type = typename container_type::mapped_type;
106 using value_type = typename container_type::value_type;
107 using iterator = typename container_type::iterator;
108 using const_iterator = typename container_type::const_iterator;
109
111 virtual ~BondedInteractionsMap() = default;
112
113 iterator begin() { return m_params.begin(); }
114 iterator end() { return m_params.end(); }
115 const_iterator begin() const { return m_params.begin(); }
116 const_iterator end() const { return m_params.end(); }
117
118 void insert(key_type const &key, mapped_type const &ptr) {
119 if (m_params.contains(key)) {
120 deactivate_bond(m_params.at(key));
121 }
122 activate_bond(ptr);
123 next_key = std::max(next_key, key + 1);
124 m_params[key] = ptr;
125 on_ia_change();
126 }
128 activate_bond(ptr);
129 auto const key = next_key++;
130 m_params[key] = ptr;
131 on_ia_change();
132 return key;
133 }
134 auto erase(key_type const &key) {
135 if (m_params.contains(key)) {
136 deactivate_bond(m_params.at(key));
137 }
138 auto &&obj = m_params.erase(key);
139 on_ia_change();
140 return obj;
141 }
142 virtual void activate_bond(mapped_type const &ptr);
143 virtual void deactivate_bond(mapped_type const &ptr);
144 mapped_type at(key_type const &key) const { return m_params.at(key); }
145 auto count(key_type const &key) const { return m_params.count(key); }
146 bool contains(key_type const &key) const { return m_params.contains(key); }
147 bool empty() const { return m_params.empty(); }
148 auto size() const { return m_params.size(); }
149 auto get_next_key() const { return next_key; }
150 auto get_zero_based_type(int bond_id) const {
151 return contains(bond_id) ? at(bond_id)->which() : 0;
152 }
153 auto get_n_thermalized_bonds() const { return n_thermalized_bonds; }
154#ifdef BOND_CONSTRAINT
155 auto get_n_rigid_bonds() const { return n_rigid_bonds; }
156#endif
157 std::optional<key_type> find_bond_id(mapped_type const &bond) const {
158 for (auto const &kv : m_params) {
159 if (kv.second == bond) {
160 return kv.first;
161 }
162 }
163 return std::nullopt;
164 }
165
166 /**
167 * @brief Calculate the maximal cutoff of bonded interactions, required to
168 * determine the cell size for communication.
169 *
170 * Bond angle and dihedral potentials do not contain a cutoff intrinsically.
171 * The cutoff for these potentials depends on the bond length potentials
172 * (it is assumed that particles participating in a bond angle or dihedral
173 * potential are bound to each other by some bond length potential). For bond
174 * angle potentials nothing has to be done. For dihedral potentials the cutoff
175 * is set to twice the maximal cutoff because the particle in which the bond
176 * is stored is only bonded to the first two partners, one of which has an
177 * additional bond to the third partner.
178 */
179 double maximal_cutoff() const;
180
181 /**
182 * @brief Checks both particles for a specific bond, even on ghost particles.
183 *
184 * @param p particle to check for the bond
185 * @param p_partner possible bond partner
186 * @tparam BondType Bond type to check for. Must be of one of the types in
187 * @ref Bonded_IA_Parameters.
188 */
189 template <typename BondType>
190 bool pair_bond_exists_on(Particle const &p, Particle const &p_partner) const {
191 auto const &bonds = p.bonds();
192 return std::any_of(
193 bonds.begin(), bonds.end(),
194 [this, partner_id = p_partner.id()](BondView const &bond) {
195 auto const &bond_ptr = at(bond.bond_id());
196 return (boost::get<BondType>(bond_ptr.get()) != nullptr) and
197 (bond.partner_ids()[0] == partner_id);
198 });
199 }
200
201 /**
202 * @brief Checks both particles for a specific bond, even on ghost particles.
203 *
204 * @param p1 particle on which the bond may be stored
205 * @param p2 particle on which the bond may be stored
206 * @tparam BondType Bond type to check for.
207 */
208 template <typename BondType>
209 bool pair_bond_exists_between(Particle const &p1, Particle const &p2) const {
210 if (&p1 == &p2)
211 return false;
212
213 // Check if particles have bonds and search for the bond of interest.
214 // Could be saved on either particle, so we need to check both.
215 return pair_bond_exists_on<BondType>(p1, p2) or
216 pair_bond_exists_on<BondType>(p2, p1);
217 }
218
219 void on_ia_change();
220
221private:
222 container_type m_params = {};
223 key_type next_key = static_cast<key_type>(0);
224 int n_thermalized_bonds = 0;
225#ifdef BOND_CONSTRAINT
226 int n_rigid_bonds = 0;
227#endif
228};
229
230/** @brief Get the number of bonded partners for the specified bond. */
231inline int number_of_partners(Bonded_IA_Parameters const &iaparams) {
232 return boost::apply_visitor(BondNumPartners(), iaparams);
233}
Common code for functions calculating angle forces.
Routines to calculate the angle energy or/and and force for a particle triple using the potential des...
Routines to calculate the angle energy or/and and force for a particle triple using the potential des...
Routines to calculate the angle energy or/and and force for a particle triple using the potential des...
Routines to calculate the bonded Coulomb potential between particle pairs.
Routines to calculate the short-range part of the bonded Coulomb potential between particle pairs.
static constexpr double BONDED_INACTIVE_CUTOFF
int number_of_partners(Bonded_IA_Parameters const &iaparams)
Get the number of bonded partners for the specified bond.
boost::variant< NoneBond, FeneBond, HarmonicBond, QuarticBond, BondedCoulomb, BondedCoulombSR, AngleHarmonicBond, AngleCosineBond, AngleCossquareBond, DihedralBond, TabulatedDistanceBond, TabulatedAngleBond, TabulatedDihedralBond, ThermalizedBond, RigidBond, IBMTriel, IBMVolCons, IBMTribend, OifGlobalForcesBond, OifLocalForcesBond, VirtualBond > Bonded_IA_Parameters
Variant in which to store the parameters of an individual bonded interaction.
Routines to calculate the energy and/or force for particle bonds, angles and dihedrals via interpolat...
Visitor to get the number of bound partners from the bond parameter variant.
int operator()(T const &) const
Immutable view on a bond.
Definition BondList.hpp:44
container for bonded interactions.
typename container_type::key_type key_type
std::optional< key_type > find_bond_id(mapped_type const &bond) const
virtual void activate_bond(mapped_type const &ptr)
auto count(key_type const &key) const
key_type insert(mapped_type const &ptr)
virtual void deactivate_bond(mapped_type const &ptr)
mapped_type at(key_type const &key) const
auto get_zero_based_type(int bond_id) const
bool contains(key_type const &key) const
virtual ~BondedInteractionsMap()=default
void insert(key_type const &key, mapped_type const &ptr)
typename container_type::const_iterator const_iterator
bool pair_bond_exists_between(Particle const &p1, Particle const &p2) const
Checks both particles for a specific bond, even on ghost particles.
bool pair_bond_exists_on(Particle const &p, Particle const &p_partner) const
Checks both particles for a specific bond, even on ghost particles.
double maximal_cutoff() const
Calculate the maximal cutoff of bonded interactions, required to determine the cell size for communic...
BondedInteractionsMap()=default
typename container_type::mapped_type mapped_type
auto erase(key_type const &key)
typename container_type::iterator iterator
typename container_type::value_type value_type
Abstract class that represents a component of the system.
Routines to calculate the dihedral energy or/and force for a particle quadruple.
Routines to calculate the FENE potential between particle pairs.
Routines to calculate the harmonic bond potential between particle pairs.
Routines to calculate the OIF local forces for a particle quadruple (two neighboring triangles with c...
Routines to calculate the quartic potential between particle pairs.
Definition of the rigid bond data type for the Rattle algorithm.
Parameters for three-body angular potential (cosine).
Parameters for three-body angular potential (cossquare).
Parameters for three-body angular potential (harmonic).
Parameters for Coulomb bond short-range Potential.
Parameters for Coulomb bond Potential.
Parameters for four-body angular potential (dihedral-angle potentials).
Definition dihedral.hpp:41
Parameters for FENE bond Potential.
Definition fene.hpp:38
Parameters for harmonic bond Potential.
Definition harmonic.hpp:37
Parameters for IBM tribend.
Parameters for IBM elastic triangle (triel)
Definition ibm_triel.hpp:35
Parameters for IBM volume conservation bond.
Interaction type for unused bonded interaction slots.
double cutoff() const
static constexpr int num
Parameters for OIF global forces.
Parameters for OIF local forces.
Struct holding all information for one particle.
Definition Particle.hpp:395
auto const & bonds() const
Definition Particle.hpp:428
auto const & id() const
Definition Particle.hpp:414
Parameters for quartic bond Potential.
Definition quartic.hpp:36
Parameters for the rigid_bond/SHAKE/RATTLE ALGORITHM.
Parameters for 3-body tabulated potential.
Parameters for 4-body tabulated potential.
Parameters for 2-body tabulated potential.
Parameters for Thermalized bond.
Interaction type for virtual bonds.
static constexpr int num
Routines to thermalize the center of mass and distance of a particle pair.