ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
ShapeBasedConstraint.cpp
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
21
22#include "BoxGeometry.hpp"
23#include "Observable_stat.hpp"
25#include "communication.hpp"
26#include "config/config.hpp"
27#include "dpd.hpp"
28#include "energy_inline.hpp"
29#include "errorhandling.hpp"
30#include "forces_inline.hpp"
32#include "system/System.hpp"
33#include "thermostat.hpp"
34
35#include <utils/Vector.hpp>
36
37#include <boost/mpi/collectives.hpp>
38
39#include <algorithm>
40#include <functional>
41#include <limits>
42#include <numeric>
43
44namespace Constraints {
45/** Check if a non-bonded interaction is defined */
46static bool is_active(IA_parameters const &data) {
47 return data.max_cut != INACTIVE_CUTOFF;
48}
49
51 part_rep.type() = type;
52 m_system.lock()->nonbonded_ias->make_particle_type_exist(type);
53}
54
55IA_parameters const &ShapeBasedConstraint::get_ia_param(int type) const {
56 return m_system.lock()->nonbonded_ias->get_ia_param(type, part_rep.type());
57}
58
60 return all_reduce(comm_cart, m_local_force, std::plus<>());
61}
62
64 return all_reduce(comm_cart, m_outer_normal_force, std::plus<double>());
65}
66
68 ParticleRange const &particles) const {
69 auto global_mindist = std::numeric_limits<double>::infinity();
70
71 auto const local_mindist = std::accumulate(
72 particles.begin(), particles.end(),
73 std::numeric_limits<double>::infinity(),
74 [this, &box_geo](double min, Particle const &p) {
75 auto const &ia_params = get_ia_param(p.type());
76 if (is_active(ia_params)) {
77 double dist;
78 Utils::Vector3d vec;
79 m_shape->calculate_dist(box_geo.folded_position(p.pos()), dist, vec);
80 return std::min(min, dist);
81 }
82 return min;
83 });
84 boost::mpi::reduce(comm_cart, local_mindist, global_mindist,
85 boost::mpi::minimum<double>(), 0);
86 return global_mindist;
87}
88
89ParticleForce ShapeBasedConstraint::force(Particle const &p,
90 Utils::Vector3d const &folded_pos,
91 double) {
92 ParticleForce pf{};
93 auto const &ia_params = get_ia_param(p.type());
94
95 if (is_active(ia_params)) {
96 double dist = 0.;
97 Utils::Vector3d dist_vec;
98 m_shape->calculate_dist(folded_pos, dist, dist_vec);
99 auto &system = *m_system.lock();
100 auto const coulomb_kernel = system.coulomb.pair_force_kernel();
101
102#ifdef DPD
103 Utils::Vector3d dpd_force{};
104#endif
105 Utils::Vector3d outer_normal_vec{};
106
107 if (dist > 0) {
108 outer_normal_vec = -dist_vec / dist;
109 pf = calc_central_radial_force(ia_params, dist_vec, dist) +
110#ifdef THOLE
111 thole_pair_force(p, part_rep, ia_params, dist_vec, dist,
112 *system.bonded_ias, get_ptr(coulomb_kernel)) +
113#endif
114 calc_non_central_force(p, part_rep, ia_params, dist_vec, dist);
115
116#ifdef DPD
117 if (system.thermostat->thermo_switch & THERMO_DPD) {
118 dpd_force = dpd_pair_force(p, part_rep, *system.thermostat->dpd,
119 *system.box_geo, ia_params, dist_vec, dist,
120 dist * dist);
121 // Additional use of DPD here requires counter increase
122 system.thermostat->dpd->rng_increment();
123 }
124#endif
125 } else if (m_penetrable && (dist <= 0)) {
126 if ((!m_only_positive) && (dist < 0)) {
127 pf = calc_central_radial_force(ia_params, dist_vec, -dist) +
128#ifdef THOLE
129 thole_pair_force(p, part_rep, ia_params, dist_vec, -dist,
130 *system.bonded_ias, get_ptr(coulomb_kernel)) +
131#endif
132 calc_non_central_force(p, part_rep, ia_params, dist_vec, -dist);
133
134#ifdef DPD
135 if (system.thermostat->thermo_switch & THERMO_DPD) {
136 dpd_force = dpd_pair_force(p, part_rep, *system.thermostat->dpd,
137 *system.box_geo, ia_params, dist_vec, dist,
138 dist * dist);
139 // Additional use of DPD here requires counter increase
140 system.thermostat->dpd->rng_increment();
141 }
142#endif
143 }
144 } else {
145 runtimeErrorMsg() << "Constraint violated by particle " << p.id()
146 << " dist " << dist;
147 }
148
149#ifdef ROTATION
150 part_rep.torque() += calc_opposing_force(pf, dist_vec).torque;
151#endif
152#ifdef DPD
153 pf.f += dpd_force;
154#endif
155 m_local_force -= pf.f;
156 m_outer_normal_force -= outer_normal_vec * pf.f;
157 }
158 return pf;
159}
160
161void ShapeBasedConstraint::add_energy(const Particle &p,
162 const Utils::Vector3d &folded_pos, double,
163 Observable_stat &obs_energy) const {
164 double energy = 0.0;
165 auto const &ia_params = get_ia_param(p.type());
166
167 if (is_active(ia_params)) {
168 auto &system = *m_system.lock();
169 auto const coulomb_kernel = system.coulomb.pair_energy_kernel();
170 double dist = 0.0;
171 Utils::Vector3d vec;
172 m_shape->calculate_dist(folded_pos, dist, vec);
173 if (dist > 0.) {
174 energy = calc_non_bonded_pair_energy(p, part_rep, ia_params, vec, dist,
175 *system.bonded_ias,
176 get_ptr(coulomb_kernel));
177 } else if (dist <= 0. and m_penetrable) {
178 if (!m_only_positive and dist < 0.) {
179 energy = calc_non_bonded_pair_energy(p, part_rep, ia_params, vec, -dist,
180 *system.bonded_ias,
181 get_ptr(coulomb_kernel));
182 }
183 } else {
184 runtimeErrorMsg() << "Constraint violated by particle " << p.id();
185 }
186 }
187 // NOLINTNEXTLINE(clang-analyzer-cplusplus.NewDeleteLeaks)
188 if (part_rep.type() >= 0) {
190 p.type(), part_rep.type(), p.mol_id(), part_rep.mol_id(), energy);
191 }
192}
193} // namespace Constraints
@ THERMO_DPD
Vector implementation and trait types for boost qvm interoperability.
Data structures for bonded interactions.
double min_dist(BoxGeometry const &box_geo, ParticleRange const &particles) const
Calculate the minimum distance between all particle pairs.
Observable for the pressure and energy.
void add_non_bonded_contribution(int type1, int type2, int molid1, int molid2, std::span< const double > data)
A range of particles.
boost::mpi::communicator comm_cart
The communicator.
This file contains the defaults for ESPResSo.
const T * get_ptr(std::optional< T > const &opt)
Utils::Vector3d dpd_pair_force(DPDParameters const &params, Utils::Vector3d const &v, double dist, Utils::Vector3d const &noise)
Definition dpd.cpp:85
Routines to use DPD as thermostat or pair force .
Energy calculation.
double calc_non_bonded_pair_energy(Particle const &p1, Particle const &p2, IA_parameters const &ia_params, Utils::Vector3d const &d, double const dist, BondedInteractionsMap const &bonded_ias, Coulomb::ShortRangeEnergyKernel::kernel_type const *coulomb_kernel)
Calculate non-bonded energies between a pair of particles.
This file contains the errorhandling code for severe errors, like a broken bond or illegal parameter ...
#define runtimeErrorMsg()
Force calculation.
ParticleForce calc_central_radial_force(IA_parameters const &ia_params, Utils::Vector3d const &d, double const dist)
ParticleForce calc_opposing_force(ParticleForce const &pf, Utils::Vector3d const &d)
ParticleForce calc_non_central_force(Particle const &p1, Particle const &p2, IA_parameters const &ia_params, Utils::Vector3d const &d, double const dist)
static bool is_active(IA_parameters const &data)
Check if a non-bonded interaction is defined.
Various procedures concerning interactions between particles.
constexpr double INACTIVE_CUTOFF
Cutoff for deactivated interactions.
Parameters for non-bonded interactions.
double max_cut
maximal cutoff for this pair of particle types.
Force information on a particle.
Definition Particle.hpp:290
Utils::Vector3d torque
torque.
Definition Particle.hpp:320
Struct holding all information for one particle.
Definition Particle.hpp:395
auto const & type() const
Definition Particle.hpp:418
auto const & mol_id() const
Definition Particle.hpp:416
auto const & id() const
Definition Particle.hpp:414
Utils::Vector3d thole_pair_force(Particle const &p1, Particle const &p2, IA_parameters const &ia_params, Utils::Vector3d const &d, double dist, BondedInteractionsMap const &bonded_ias, Coulomb::ShortRangeForceKernel::kernel_type const *kernel)
Calculate Thole force.
Definition thole.hpp:44