Skip to content
Permalink
Browse files

Merge pull request #2542 from fweik/bibo

core: Removed duplication of pair bond force calc.
  • Loading branch information...
fweik committed Mar 14, 2019
2 parents 507a13a + d1142ca commit 53d18aa98fa126e2a36ee0c3af638df82378c4dd
Showing with 1,083 additions and 1,276 deletions.
  1. +9 −8 doc/sphinx/inter_bonded.rst
  2. +19 −0 src/core/TabulatedPotential.hpp
  3. +120 −0 src/core/bonded_interactions/angle_common.hpp
  4. +66 −154 src/core/bonded_interactions/angle_cosine.hpp
  5. +55 −139 src/core/bonded_interactions/angle_cossquare.hpp
  6. +62 −164 src/core/bonded_interactions/angle_harmonic.hpp
  7. +23 −16 src/core/bonded_interactions/bonded_coulomb.hpp
  8. +24 −16 src/core/bonded_interactions/bonded_coulomb_p3m_sr.hpp
  9. +58 −42 src/core/bonded_interactions/bonded_interaction_data.hpp
  10. +183 −205 src/core/bonded_interactions/bonded_tab.hpp
  11. +11 −9 src/core/bonded_interactions/dihedral.hpp
  12. +24 −16 src/core/bonded_interactions/fene.hpp
  13. +24 −16 src/core/bonded_interactions/harmonic.hpp
  14. +23 −16 src/core/bonded_interactions/harmonic_dumbbell.hpp
  15. +22 −16 src/core/bonded_interactions/quartic.hpp
  16. +11 −12 src/core/bonded_interactions/thermalized_bond.hpp
  17. +22 −8 src/core/bonded_interactions/umbrella.hpp
  18. +1 −1 src/core/electrostatics_magnetostatics/p3m-dipolar.hpp
  19. +1 −1 src/core/electrostatics_magnetostatics/p3m.hpp
  20. +69 −79 src/core/forces_inline.hpp
  21. +0 −1 src/core/nonbonded_interactions/nonbonded_tab.hpp
  22. +31 −217 src/core/pressure_inline.hpp
  23. +4 −4 src/core/utils.hpp
  24. +7 −0 src/core/utils/linear_interpolation.hpp
  25. +2 −2 src/python/espressomd/integrate.pyx
  26. +0 −1 src/python/espressomd/interactions.pxd
  27. +16 −1 src/python/espressomd/interactions.pyx
  28. +1 −0 testsuite/python/CMakeLists.txt
  29. +79 −47 testsuite/python/interactions_bond_angle.py
  30. +67 −82 testsuite/python/interactions_bonded.py
  31. +49 −3 testsuite/python/interactions_dihedral.py
@@ -93,7 +93,7 @@ Harmonic Dumbbell Bond
Requires ``ROTATION`` feature.


A harmonic bond can be instantiated via
A harmonic Dumbbell bond can be instantiated via
:class:`espressomd.interactions.HarmonicDumbbellBond`::

from espressomd.interactions import HarmonicDumbbellBond
@@ -104,11 +104,11 @@ This bond is similar to the normal harmonic bond in such a way that it
sets up a harmonic potential, i.e. a spring, between the two particles.
Additionally the orientation of the first particle in the bond will be aligned along
the distance vector between both particles. This alignment can be
controlled by the second harmonic constant :math:`k2`. Keep in mind that orientation will
controlled by the second harmonic constant :math:`k_2`. Keep in mind that orientation will
oscillate around the distance vector and some kind of
friction needs to be present for it to relax.

The roles of the parameters :math:`k1, r_0, r_\mathrm{cut}` are exactly the same as for the
The roles of the parameters :math:`k_1, r_0, r_\mathrm{cut}` are exactly the same as for the
harmonic bond.

..
@@ -154,7 +154,7 @@ It is given by

.. math:: V(r) = \frac{\alpha q_1 q_2}{r},

where :math:`q1` and :math:`q2` are the charges of the bound particles and :math:`\alpha` is the
where :math:`q_1` and :math:`q_2` are the charges of the bound particles and :math:`\alpha` is the
Coulomb prefactor. This interaction has no cutoff and acts independently of other
Coulomb interactions.

@@ -575,6 +575,7 @@ The parameter ``bond_angle`` is a bond type identifier of three possible bond-an


:class:`espressomd.interactions.AngleHarmonic`

A classical harmonic potential of the form:

.. math:: V(\phi) = \frac{K}{2} \left(\phi - \phi_0\right)^2.
@@ -590,7 +591,7 @@ The parameter ``bond_angle`` is a bond type identifier of three possible bond-an
:math:`\phi=\phi_0+\pi` and accordingly a discontinuity in the
force, and should therefore be used with caution.

example ::
Example::

>>> angle_harmonic = AngleHarmonic(bend=1.0, phi0=np.pi)
>>> system.bonded_inter.add(angle_harmonic)
@@ -602,7 +603,7 @@ The parameter ``bond_angle`` is a bond type identifier of three possible bond-an

Cosine bond angle potential of the form:

.. math:: V(\phi) = K \left[1 - \cos(\phi - \phi0)\right]
.. math:: V(\phi) = K \left[1 - \cos(\phi - \phi_0)\right]

:math:`K` is the bending constant,
and the optional parameter :math:`\phi_0` is the equilibrium bond angle in
@@ -615,7 +616,7 @@ The parameter ``bond_angle`` is a bond type identifier of three possible bond-an
(both are :math:`1/2(\phi-\phi_0)^2` in leading order), but it is
periodic and smooth for all angles :math:`\phi`.

example ::
Example::

>>> angle_cosine = AngleCosine(bend=1.0, phi0=np.pi)
>>> system.bonded_inter.add(angle_cosine)
@@ -631,7 +632,7 @@ The parameter ``bond_angle`` is a bond type identifier of three possible bond-an
potential is :math:`1/8(\phi-\phi_0)^4` around :math:`\phi_0`, and
therefore much flatter than the two potentials before.

example ::
Example::

>>> angle_cossquare = AngleCossquare(bend=1.0, phi0=np.pi)
>>> system.bonded_inter.add(angle_cossquare)
@@ -29,19 +29,38 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
#include <cassert>
#include <vector>

/** Evaluate forces and energies using a custom potential profile.
*
* Forces and energies are evaluated by linear interpolation. The curves
* @ref force_tab and @ref energy_tab must be sampled uniformly between
* @ref minval and @ref maxval.
*/
struct TabulatedPotential {
/** Position on the x-axis of the first tabulated value. */
double minval = -1.0;
/** Position on the x-axis of the last tabulated value. */
double maxval = -1.0;
/** %Distance on the x-axis between tabulated values. */
double invstepsize = 0.0;
/** Tabulated forces. */
std::vector<double> force_tab;
/** Tabulated energies. */
std::vector<double> energy_tab;

/** Evaluate the force at position @p x.
* @param x Bond length/angle
* @return Interpolated force.
*/
double force(double x) const {
using boost::algorithm::clamp;
return Utils::linear_interpolation(force_tab, invstepsize, minval,
clamp(x, minval, maxval));
}

/** Evaluate the energy at position @p x.
* @param x Bond length/angle
* @return Interpolated energy.
*/
double energy(double x) const {
using boost::algorithm::clamp;
return Utils::linear_interpolation(energy_tab, invstepsize, minval,
@@ -0,0 +1,120 @@
/*
Copyright (C) 2010-2019 The ESPResSo project
Copyright (C) 2002-2010
Max-Planck-Institute for Polymer Research, Theory Group
This file is part of ESPResSo.
ESPResSo is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
ESPResSo is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef ANGLE_COMMON_H
#define ANGLE_COMMON_H
/** \file
* Common code for functions calculating angle forces.
*/

#include "grid.hpp"
#include <tuple>

inline std::tuple<Vector3d, Vector3d, double, double, double>
calc_vectors_and_cosine(Vector3d const &r_mid, Vector3d const &r_left,
Vector3d const &r_right, bool sanitize_cosine = false) {
/* vector from p_left to p_mid */
auto vec1 = get_mi_vector(r_mid, r_left);
auto const d1i = 1.0 / vec1.norm();
vec1 *= d1i;
/* vector from p_mid to p_right */
auto vec2 = get_mi_vector(r_right, r_mid);
auto const d2i = 1.0 / vec2.norm();
vec2 *= d2i;
/* cosine of the angle between vec1 and vec2 */
auto cosine = vec1 * vec2;
if (sanitize_cosine) {
if (cosine > TINY_COS_VALUE)
cosine = TINY_COS_VALUE;
if (cosine < -TINY_COS_VALUE)
cosine = -TINY_COS_VALUE;
}
return std::make_tuple(vec1, vec2, d1i, d2i, cosine);
}

/** Compute a three-body angle interaction force.
* @param[in] r_mid Position of second/middle particle.
* @param[in] r_left Position of first/left particle.
* @param[in] r_right Position of third/right particle.
* @param[in] forceFactor Angle bending constant.
* @param[in] sanitize_cosine Sanitize the cosine of the angle.
* @param[out] force1 Force on particle 1.
* @param[out] force2 Force on particle 2.
*/
template <typename ForceFactor>
void calc_angle_generic_force(Vector3d const &r_mid, Vector3d const &r_left,
Vector3d const &r_right, ForceFactor forceFactor,
double force1[3], double force2[3],
bool sanitize_cosine) {
Vector3d vec1, vec2;
double d1i, d2i, cosine;
std::tie(vec1, vec2, d1i, d2i, cosine) =
calc_vectors_and_cosine(r_mid, r_left, r_right, sanitize_cosine);
/* force factor */
auto const fac = forceFactor(cosine);
/* force calculation */
auto const f1 = fac * (cosine * vec1 - vec2) * d1i;
auto const f2 = fac * (cosine * vec2 - vec1) * d2i;
for (int j = 0; j < 3; j++) {
force1[j] = (f1 - f2)[j];
force2[j] = -f1[j];
}
}

/** Compute the forces of a three-body bonded potential.
* @param[in] r_mid Position of second/middle particle.
* @param[in] r_left Position of first/left particle.
* @param[in] r_right Position of third/right particle.
* @param[in] forceFactor Angle bending constant.
* @param[in] sanitize_cosine Sanitize the cosine of the angle.
* @return Forces on particles 1, 2 and 3.
*/
template <typename ForceFactor>
std::tuple<Vector3d, Vector3d, Vector3d>
calc_angle_generic_3body_forces(Vector3d const &r_mid, Vector3d const &r_left,
Vector3d const &r_right,
ForceFactor forceFactor, bool sanitize_cosine) {
auto const vec21 = get_mi_vector(r_left, r_mid);
auto const vec31 = get_mi_vector(r_right, r_mid);
auto const vec21_sqr = vec21.norm2();
auto const vec31_sqr = vec31.norm2();
auto const vec21_len = sqrt(vec21_sqr);
auto const vec31_len = sqrt(vec31_sqr);
auto cos_phi = (vec21 * vec31) / (vec21_len * vec31_len);
if (sanitize_cosine) {
if (cos_phi > TINY_COS_VALUE)
cos_phi = TINY_COS_VALUE;
if (cos_phi < -TINY_COS_VALUE)
cos_phi = -TINY_COS_VALUE;
}
auto sin_phi = sqrt(1.0 - Utils::sqr(cos_phi));
/* force factor */
auto const fac = forceFactor(cos_phi, sin_phi);
/* force calculation */
auto const fj = vec31 / (vec21_len * vec31_len) - cos_phi * vec21 / vec21_sqr;
auto const fk = vec21 / (vec21_len * vec31_len) - cos_phi * vec31 / vec31_sqr;
// note that F1 = -(F2 + F3) in analytical case
auto force1 = -fac * (fj + fk);
auto force2 = fac * fj;
auto force3 = fac * fk;
return std::make_tuple(force1, force2, force3);
}

#endif /* ANGLE_COMMON_H */
Oops, something went wrong.

0 comments on commit 53d18aa

Please sign in to comment.
You can’t perform that action at this time.
You signed in with another tab or window. Reload to refresh your session. You signed out in another tab or window. Reload to refresh your session.