ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
bonded_tab.hpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-2022 The ESPResSo project
3 * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
4 * Max-Planck-Institute for Polymer Research, Theory Group
5 *
6 * This file is part of ESPResSo.
7 *
8 * ESPResSo is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * ESPResSo is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with this program. If not, see <http://www.gnu.org/licenses/>.
20 */
21
22#pragma once
23
24/** \file
25 * Routines to calculate the energy and/or force for particle bonds, angles
26 * and dihedrals via interpolation of lookup tables.
27 */
28
29#include "config/config.hpp"
30
32#include "angle_common.hpp"
34
35#include <utils/Vector.hpp>
36#include <utils/math/sqr.hpp>
37
38#include <cassert>
39#include <cmath>
40#include <memory>
41#include <numbers>
42#include <optional>
43#include <tuple>
44#include <vector>
45
46/** Base class for n-body tabulated potential (n=2,3,4). */
48 std::shared_ptr<TabulatedPotential> pot;
49
50 /** Set the parameters of a bonded tabulated potential.
51 * ia_params and force/energy tables are communicated to each node.
52 *
53 * @param min @copybrief TabulatedPotential::minval
54 * @param max @copybrief TabulatedPotential::maxval
55 * @param energy @copybrief TabulatedPotential::energy_tab
56 * @param force @copybrief TabulatedPotential::force_tab
57 */
58 TabulatedBond(double min, double max, std::vector<double> const &energy,
59 std::vector<double> const &force) {
60 pot = std::make_shared<TabulatedPotential>(min, max, force, energy);
61 }
62};
63
64/** Parameters for 2-body tabulated potential. */
66 double cutoff() const { return assert(pot), pot->cutoff(); }
67
68 static constexpr int num = 1;
69
70 TabulatedDistanceBond(double min, double max,
71 std::vector<double> const &energy,
72 std::vector<double> const &force)
73 : TabulatedBond(min, max, energy, force) {
74 this->pot->minval = min;
75 this->pot->maxval = max;
76 }
77
78 std::optional<Utils::Vector3d> force(Utils::Vector3d const &dx) const;
79 std::optional<double> energy(Utils::Vector3d const &dx) const;
80};
81
82/** Parameters for 3-body tabulated potential. */
84 double cutoff() const { return 0.; }
85
86 static constexpr int num = 2;
87
88 TabulatedAngleBond(double min, double max, std::vector<double> const &energy,
89 std::vector<double> const &force)
90 : TabulatedBond(min, max, energy, force) {
91 this->pot->minval = 0.;
92 this->pot->maxval = std::numbers::pi + ROUND_ERROR_PREC;
93 }
94
95 std::tuple<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>
96 forces(Utils::Vector3d const &vec1, Utils::Vector3d const &vec2) const;
97 double energy(Utils::Vector3d const &vec1, Utils::Vector3d const &vec2) const;
98};
99
100/** Parameters for 4-body tabulated potential. */
102 double cutoff() const { return 0.; }
103
104 static constexpr int num = 3;
105
106 TabulatedDihedralBond(double min, double max,
107 std::vector<double> const &energy,
108 std::vector<double> const &force)
109 : TabulatedBond(min, max, energy, force) {
110 this->pot->minval = 0.;
111 this->pot->maxval = 2. * std::numbers::pi + ROUND_ERROR_PREC;
112 }
113
114 std::optional<std::tuple<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d,
116 forces(Utils::Vector3d const &v12, Utils::Vector3d const &v23,
117 Utils::Vector3d const &v34) const;
118 std::optional<double> energy(Utils::Vector3d const &v12,
119 Utils::Vector3d const &v23,
120 Utils::Vector3d const &v34) const;
121};
122
123/** Compute a tabulated bond length force.
124 *
125 * The force acts in the direction of the connecting vector between the
126 * particles. For distances smaller than the tabulated range it uses a linear
127 * extrapolation based on the first two tabulated force values.
128 *
129 * @param[in] dx Distance between the particles.
130 */
131inline std::optional<Utils::Vector3d>
133 auto const dist = dx.norm();
134
135 if (dist < pot->cutoff()) {
136 auto const fac = pot->force(dist) / dist;
137 return fac * dx;
138 }
139 return {};
140}
141
142/** Compute a tabulated bond length energy.
143 *
144 * For distances smaller than the tabulated range it uses a quadratic
145 * extrapolation based on the first two tabulated force values and the first
146 * tabulated energy value.
147 *
148 * @param[in] dx Distance between the particles.
149 */
150inline std::optional<double>
152 auto const dist = dx.norm();
153
154 if (dist < pot->cutoff()) {
155 return pot->energy(dist);
156 }
157 return {};
158}
159
160/** Compute the three-body angle interaction force.
161 * @param[in] vec1 Vector from central particle to left particle.
162 * @param[in] vec2 Vector from central particle to right particle.
163 * @return Forces on the second, first and third particles, in that order.
164 */
165inline std::tuple<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>
167 Utils::Vector3d const &vec2) const {
168
169 auto forceFactor = [this](double const cos_phi) {
170 auto const sin_phi = sqrt(1 - Utils::sqr(cos_phi));
171#ifdef TABANGLEMINUS
172 auto const phi = acos(-cos_phi);
173#else
174 auto const phi = acos(cos_phi);
175#endif
176 auto const tab_pot = pot;
177 auto const gradient = tab_pot->force(phi);
178 return -gradient / sin_phi;
179 };
180
181 return angle_generic_force(vec1, vec2, forceFactor, true);
182}
183
184/** Compute the three-body angle interaction energy.
185 * It is assumed that the potential is tabulated
186 * for all angles between 0 and Pi.
187 *
188 * @param[in] vec1 Vector from central particle to left particle.
189 * @param[in] vec2 Vector from central particle to right particle.
190 */
192 Utils::Vector3d const &vec2) const {
193 auto const cos_phi = calc_cosine(vec1, vec2, true);
194 /* calculate phi */
195#ifdef TABANGLEMINUS
196 auto const phi = acos(-cos_phi);
197#else
198 auto const phi = acos(cos_phi);
199#endif
200 return pot->energy(phi);
201}
202
203/** Compute the four-body dihedral interaction force.
204 * The forces have a singularity at @f$ \phi = 0 @f$ and @f$ \phi = \pi @f$
205 * (see @cite swope92a page 592).
206 *
207 * @param[in] v12 Vector from @p p1 to @p p2
208 * @param[in] v23 Vector from @p p2 to @p p3
209 * @param[in] v34 Vector from @p p3 to @p p4
210 * @return the forces on @p p2, @p p1, @p p3
211 */
212inline std::optional<std::tuple<Utils::Vector3d, Utils::Vector3d,
215 Utils::Vector3d const &v23,
216 Utils::Vector3d const &v34) const {
217 /* vectors for dihedral angle calculation */
218 Utils::Vector3d v12Xv23, v23Xv34;
219 double l_v12Xv23, l_v23Xv34;
220 /* dihedral angle, cosine of the dihedral angle */
221 double phi, cos_phi;
222
223 /* dihedral angle */
224 auto const angle_is_undefined = calc_dihedral_angle(
225 v12, v23, v34, v12Xv23, l_v12Xv23, v23Xv34, l_v23Xv34, cos_phi, phi);
226 /* dihedral angle not defined - force zero */
227 if (angle_is_undefined) {
228 return {};
229 }
230
231 auto const f1 = (v23Xv34 - cos_phi * v12Xv23) / l_v12Xv23;
232 auto const f4 = (v12Xv23 - cos_phi * v23Xv34) / l_v23Xv34;
233
234 auto const v23Xf1 = vector_product(v23, f1);
235 auto const v23Xf4 = vector_product(v23, f4);
236 auto const v34Xf4 = vector_product(v34, f4);
237 auto const v12Xf1 = vector_product(v12, f1);
238
239 /* table lookup */
240 auto const fac = pot->force(phi);
241
242 /* store dihedral forces */
243 auto const force1 = fac * v23Xf1;
244 auto const force2 = fac * (v34Xf4 - v12Xf1 - v23Xf1);
245 auto const force3 = fac * (v12Xf1 - v23Xf4 - v34Xf4);
246
247 return std::make_tuple(force2, force1, force3, -(force2 + force1 + force3));
248}
249
250/** Compute the four-body dihedral interaction energy.
251 * The energy doesn't have any singularity if the angle phi is well-defined.
252 *
253 * @param[in] v12 Vector from @p p1 to @p p2
254 * @param[in] v23 Vector from @p p2 to @p p3
255 * @param[in] v34 Vector from @p p3 to @p p4
256 */
257inline std::optional<double>
259 Utils::Vector3d const &v23,
260 Utils::Vector3d const &v34) const {
261 /* vectors for dihedral calculations. */
262 Utils::Vector3d v12Xv23, v23Xv34;
263 double l_v12Xv23, l_v23Xv34;
264 /* dihedral angle, cosine of the dihedral angle */
265 double phi, cos_phi;
266 auto const angle_is_undefined = calc_dihedral_angle(
267 v12, v23, v34, v12Xv23, l_v12Xv23, v23Xv34, l_v23Xv34, cos_phi, phi);
268 /* dihedral angle not defined - energy zero */
269 if (angle_is_undefined) {
270 return {};
271 }
272
273 return pot->energy(phi);
274}
Vector implementation and trait types for boost qvm interoperability.
Common code for functions calculating angle forces.
std::tuple< Utils::Vector3d, Utils::Vector3d, Utils::Vector3d > angle_generic_force(Utils::Vector3d const &vec1, Utils::Vector3d const &vec2, ForceFactor forceFactor, bool sanitize_cosine)
Compute a three-body angle interaction force.
double calc_cosine(Utils::Vector3d const &vec1, Utils::Vector3d const &vec2, bool sanitize_cosine=false)
Compute the cosine of the angle between three particles.
T norm() const
Definition Vector.hpp:138
This file contains the defaults for ESPResSo.
#define ROUND_ERROR_PREC
Precision for capture of round off errors.
Definition config.hpp:66
Routines to calculate the dihedral energy or/and force for a particle quadruple.
bool calc_dihedral_angle(Utils::Vector3d const &a, Utils::Vector3d const &b, Utils::Vector3d const &c, Utils::Vector3d &aXb, double &l_aXb, Utils::Vector3d &bXc, double &l_bXc, double &cosphi, double &phi)
Calculates the dihedral angle between particle quadruple p1, p2, p3 and p4.
Definition dihedral.hpp:87
__device__ void vector_product(float const *a, float const *b, float *out)
VectorXd< 3 > Vector3d
Definition Vector.hpp:164
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
Definition sqr.hpp:28
Parameters for 3-body tabulated potential.
TabulatedAngleBond(double min, double max, std::vector< double > const &energy, std::vector< double > const &force)
double energy(Utils::Vector3d const &vec1, Utils::Vector3d const &vec2) const
Compute the three-body angle interaction energy.
static constexpr int num
std::tuple< Utils::Vector3d, Utils::Vector3d, Utils::Vector3d > forces(Utils::Vector3d const &vec1, Utils::Vector3d const &vec2) const
Compute the three-body angle interaction force.
double cutoff() const
Base class for n-body tabulated potential (n=2,3,4).
std::shared_ptr< TabulatedPotential > pot
TabulatedBond(double min, double max, std::vector< double > const &energy, std::vector< double > const &force)
Set the parameters of a bonded tabulated potential.
Parameters for 4-body tabulated potential.
std::optional< std::tuple< Utils::Vector3d, Utils::Vector3d, Utils::Vector3d, Utils::Vector3d > > forces(Utils::Vector3d const &v12, Utils::Vector3d const &v23, Utils::Vector3d const &v34) const
Compute the four-body dihedral interaction force.
TabulatedDihedralBond(double min, double max, std::vector< double > const &energy, std::vector< double > const &force)
static constexpr int num
double cutoff() const
std::optional< double > energy(Utils::Vector3d const &v12, Utils::Vector3d const &v23, Utils::Vector3d const &v34) const
Compute the four-body dihedral interaction energy.
Parameters for 2-body tabulated potential.
double cutoff() const
std::optional< Utils::Vector3d > force(Utils::Vector3d const &dx) const
Compute a tabulated bond length force.
static constexpr int num
std::optional< double > energy(Utils::Vector3d const &dx) const
Compute a tabulated bond length energy.
TabulatedDistanceBond(double min, double max, std::vector< double > const &energy, std::vector< double > const &force)