39#include <boost/optional.hpp>
40#include <boost/serialization/shared_ptr.hpp>
50 std::shared_ptr<TabulatedPotential>
pot;
60 TabulatedBond(
double min,
double max, std::vector<double>
const &energy,
61 std::vector<double>
const &
force) {
62 pot = std::make_shared<TabulatedPotential>(min, max,
force, energy);
66 friend boost::serialization::access;
67 template <
typename Archive>
68 void serialize(Archive &ar,
long int ) {
77 static constexpr int num = 1;
80 std::vector<double>
const &
energy,
81 std::vector<double>
const &
force)
83 this->
pot->minval = min;
84 this->
pot->maxval = max;
95 static constexpr int num = 2;
98 std::vector<double>
const &
force)
100 this->
pot->minval = 0.;
104 std::tuple<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>
113 static constexpr int num = 3;
116 std::vector<double>
const &
energy,
117 std::vector<double>
const &
force)
119 this->
pot->minval = 0.;
140inline boost::optional<Utils::Vector3d>
142 auto const dist = dx.
norm();
144 if (dist < pot->
cutoff()) {
145 auto const fac =
pot->force(dist) / dist;
159inline boost::optional<double>
161 auto const dist = dx.
norm();
163 if (dist < pot->
cutoff()) {
164 return pot->energy(dist);
174inline std::tuple<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>
178 auto forceFactor = [
this](
double const cos_phi) {
179 auto const sin_phi = sqrt(1 -
Utils::sqr(cos_phi));
181 auto const phi = acos(-cos_phi);
183 auto const phi = acos(cos_phi);
185 auto const tab_pot =
pot;
186 auto const gradient = tab_pot->force(phi);
187 return -gradient / sin_phi;
202 auto const cos_phi =
calc_cosine(vec1, vec2,
true);
205 auto const phi = acos(-cos_phi);
207 auto const phi = acos(cos_phi);
209 return pot->energy(phi);
228 double l_v12Xv23, l_v23Xv34;
234 v12, v23, v34, v12Xv23, l_v12Xv23, v23Xv34, l_v23Xv34, cos_phi, phi);
236 if (angle_is_undefined) {
240 auto const f1 = (v23Xv34 - cos_phi * v12Xv23) / l_v12Xv23;
241 auto const f4 = (v12Xv23 - cos_phi * v23Xv34) / l_v23Xv34;
249 auto const fac =
pot->force(phi);
252 auto const force1 = fac * v23Xf1;
253 auto const force2 = fac * (v34Xf4 - v12Xf1 - v23Xf1);
254 auto const force3 = fac * (v12Xf1 - v23Xf4 - v34Xf4);
256 return std::make_tuple(force2, force1, force3, -(force2 + force1 + force3));
266inline boost::optional<double>
272 double l_v12Xv23, l_v23Xv34;
276 v12, v23, v34, v12Xv23, l_v12Xv23, v23Xv34, l_v23Xv34, cos_phi, phi);
278 if (angle_is_undefined) {
282 return pot->energy(phi);
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.
This file contains the defaults for ESPResSo.
#define ROUND_ERROR_PREC
Precision for capture of round off errors.
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.
__device__ void vector_product(float const *a, float const *b, float *out)
DEVICE_QUALIFIER constexpr T pi()
Ratio of diameter and circumference of a circle.
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
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.
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.
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.
boost::optional< double > energy(Utils::Vector3d const &v12, Utils::Vector3d const &v23, Utils::Vector3d const &v34) const
Compute the four-body dihedral interaction energy.
TabulatedDihedralBond(double min, double max, std::vector< double > const &energy, std::vector< double > const &force)
boost::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.
Parameters for 2-body tabulated potential.
boost::optional< Utils::Vector3d > force(Utils::Vector3d const &dx) const
Compute a tabulated bond length force.
boost::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)