ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
dihedral.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 dihedral energy or/and
26 * force for a particle quadruple. Note that usage of dihedrals
27 * increases the interaction range of bonded interactions to 2 times
28 * the maximal bond length!
29 */
30
31#include "config/config.hpp"
32
33#include <utils/Vector.hpp>
34
35#include <cmath>
36#include <numbers>
37#include <optional>
38#include <tuple>
39
40/** Parameters for four-body angular potential (dihedral-angle potentials). */
42 double mult;
43 double bend;
44 double phase;
45
46 double cutoff() const { return 0.; }
47
48 static constexpr int num = 3;
49
50 DihedralBond(int mult, double bend, double phase) {
51 this->mult = mult;
52 this->bend = bend;
53 this->phase = phase;
54 }
55
56 std::optional<std::tuple<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d,
58 forces(Utils::Vector3d const &v12, Utils::Vector3d const &v23,
59 Utils::Vector3d const &v34) const;
60
61 std::optional<double> energy(Utils::Vector3d const &v12,
62 Utils::Vector3d const &v23,
63 Utils::Vector3d const &v34) const;
64};
65
66/**
67 * @brief Calculates the dihedral angle between particle quadruple p1, p2, p3
68 * and p4.
69 *
70 * The dihedral angle is the angle between the planes
71 * specified by the particle triples (p1,p2,p3) and (p2,p3,p4).
72 * Vectors a, b and c are the bond vectors between consecutive particles.
73 * If the a,b or b,c are parallel the dihedral angle is not defined in which
74 * case the function returns true. Calling functions should check for that.
75 *
76 * @param[in] a Vector from @p p1 to @p p2
77 * @param[in] b Vector from @p p2 to @p p3
78 * @param[in] c Vector from @p p3 to @p p4
79 * @param[out] aXb Vector product of a and b
80 * @param[out] l_aXb |aXB|
81 * @param[out] bXc Vector product of b and c
82 * @param[out] l_bXc |bXc|
83 * @param[out] cosphi Cosine of the dihedral angle
84 * @param[out] phi Dihedral angle in the range [0, pi]
85 * @return Whether the angle is undefined.
86 */
88 Utils::Vector3d const &b,
89 Utils::Vector3d const &c, Utils::Vector3d &aXb,
90 double &l_aXb, Utils::Vector3d &bXc,
91 double &l_bXc, double &cosphi, double &phi) {
92
93 /* calculate vector product a X b and b X c */
94 aXb = vector_product(a, b);
95 bXc = vector_product(b, c);
96
97 /* calculate the unit vectors */
98 l_aXb = aXb.norm();
99 l_bXc = bXc.norm();
100
101 /* catch case of undefined dihedral angle */
102 if (l_aXb <= TINY_LENGTH_VALUE || l_bXc <= TINY_LENGTH_VALUE) {
103 phi = -1.;
104 cosphi = 0.;
105 return true;
106 }
107
108 aXb /= l_aXb;
109 bXc /= l_bXc;
110
111 cosphi = aXb * bXc;
112
113 if (fabs(fabs(cosphi) - 1.) < TINY_SIN_VALUE)
114 cosphi = std::round(cosphi);
115
116 /* Calculate dihedral angle */
117 phi = acos(cosphi);
118 if ((aXb * c) < 0.)
119 phi = 2. * std::numbers::pi - phi;
120 return false;
121}
122
123/** Compute the four-body dihedral interaction force.
124 * The forces have a singularity at @f$ \phi = 0 @f$ and @f$ \phi = \pi @f$
125 * (see @cite swope92a page 592).
126 *
127 * @param[in] v12 Vector from @p p1 to @p p2
128 * @param[in] v23 Vector from @p p2 to @p p3
129 * @param[in] v34 Vector from @p p3 to @p p4
130 * @return the forces on @p p2, @p p1, @p p3
131 */
132inline std::optional<std::tuple<Utils::Vector3d, Utils::Vector3d,
135 Utils::Vector3d const &v34) const {
136 /* vectors for dihedral angle calculation */
137 Utils::Vector3d v12Xv23, v23Xv34;
138 double l_v12Xv23, l_v23Xv34;
139 /* dihedral angle, cosine of the dihedral angle */
140 double phi, cos_phi, sin_mphi_over_sin_phi;
141
142 /* dihedral angle */
143 auto const angle_is_undefined = calc_dihedral_angle(
144 v12, v23, v34, v12Xv23, l_v12Xv23, v23Xv34, l_v23Xv34, cos_phi, phi);
145 /* dihedral angle not defined - force zero */
146 if (angle_is_undefined) {
147 return {};
148 }
149
150 auto const f1 = (v23Xv34 - cos_phi * v12Xv23) / l_v12Xv23;
151 auto const f4 = (v12Xv23 - cos_phi * v23Xv34) / l_v23Xv34;
152
153 auto const v23Xf1 = vector_product(v23, f1);
154 auto const v23Xf4 = vector_product(v23, f4);
155 auto const v34Xf4 = vector_product(v34, f4);
156 auto const v12Xf1 = vector_product(v12, f1);
157
158 /* calculate force magnitude */
159 auto fac = -bend * mult;
160
161 if (fabs(sin(phi)) < TINY_SIN_VALUE) {
162 /* comes from taking the first term of the MacLaurin expansion of
163 * sin(n * phi - phi0) and sin(phi) and then making the division */
164 sin_mphi_over_sin_phi = mult * cos(mult * phi - phase) / cos_phi;
165 } else {
166 sin_mphi_over_sin_phi = sin(mult * phi - phase) / sin(phi);
167 }
168
169 fac *= sin_mphi_over_sin_phi;
170
171 /* store dihedral forces */
172 auto const force1 = fac * v23Xf1;
173 auto const force2 = fac * (v34Xf4 - v12Xf1 - v23Xf1);
174 auto const force3 = fac * (v12Xf1 - v23Xf4 - v34Xf4);
175
176 return std::make_tuple(force2, force1, force3, -(force1 + force2 + force3));
177}
178
179/** Compute the four-body dihedral interaction energy.
180 * The energy doesn't have any singularity if the angle phi is well-defined.
181 *
182 * @param[in] v12 Vector from @p p1 to @p p2
183 * @param[in] v23 Vector from @p p2 to @p p3
184 * @param[in] v34 Vector from @p p3 to @p p4
185 */
186inline std::optional<double>
188 Utils::Vector3d const &v34) const {
189 /* vectors for dihedral calculations. */
190 Utils::Vector3d v12Xv23, v23Xv34;
191 double l_v12Xv23, l_v23Xv34;
192 /* dihedral angle, cosine of the dihedral angle */
193 double phi, cos_phi;
194
195 auto const angle_is_undefined = calc_dihedral_angle(
196 v12, v23, v34, v12Xv23, l_v12Xv23, v23Xv34, l_v23Xv34, cos_phi, phi);
197 /* dihedral angle not defined - energy zero */
198 if (angle_is_undefined) {
199 return {};
200 }
201
202 return bend * (1. - cos(mult * phi - phase));
203}
Vector implementation and trait types for boost qvm interoperability.
T norm() const
Definition Vector.hpp:138
This file contains the defaults for ESPResSo.
#define TINY_LENGTH_VALUE
Tiny length cutoff.
Definition config.hpp:79
#define TINY_SIN_VALUE
Tiny angle cutoff for sinus calculations.
Definition config.hpp:71
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
Parameters for four-body angular potential (dihedral-angle potentials).
Definition dihedral.hpp:41
double phase
Definition dihedral.hpp:44
double cutoff() const
Definition dihedral.hpp:46
std::optional< double > energy(Utils::Vector3d const &v12, Utils::Vector3d const &v23, Utils::Vector3d const &v34) const
Compute the four-body dihedral interaction energy.
Definition dihedral.hpp:187
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.
Definition dihedral.hpp:134
DihedralBond(int mult, double bend, double phase)
Definition dihedral.hpp:50
double mult
Definition dihedral.hpp:42
double bend
Definition dihedral.hpp:43
static constexpr int num
Definition dihedral.hpp:48