From e945b9a27fd0ac1d556bc47a22650b6d64c2ab62 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 20 Sep 2023 17:55:36 +0200 Subject: [PATCH] core: Move BoxGeometry out of bonded IA kernels --- src/core/bonded_interactions/angle_common.hpp | 76 +++++++------- src/core/bonded_interactions/angle_cosine.hpp | 38 +++---- .../bonded_interactions/angle_cossquare.hpp | 38 +++---- .../bonded_interactions/angle_harmonic.hpp | 38 +++---- .../bonded_interactions.dox | 4 +- src/core/bonded_interactions/bonded_tab.hpp | 98 +++++++------------ src/core/bonded_interactions/dihedral.hpp | 77 ++++++--------- src/core/energy_inline.hpp | 18 ++-- src/core/forces_inline.hpp | 20 ++-- src/core/immersed_boundary/ibm_triel.cpp | 11 +-- src/core/immersed_boundary/ibm_triel.hpp | 5 +- 11 files changed, 167 insertions(+), 256 deletions(-) diff --git a/src/core/bonded_interactions/angle_common.hpp b/src/core/bonded_interactions/angle_common.hpp index 12e5726ccf5..8c12849be3f 100644 --- a/src/core/bonded_interactions/angle_common.hpp +++ b/src/core/bonded_interactions/angle_common.hpp @@ -24,52 +24,42 @@ * Common code for functions calculating angle forces. */ -#include "BoxGeometry.hpp" #include "config/config.hpp" #include +#include #include +namespace detail { +inline double sanitize_cosine(double cosine) { + if (cosine > TINY_COS_VALUE) + cosine = TINY_COS_VALUE; + if (cosine < -TINY_COS_VALUE) + cosine = -TINY_COS_VALUE; + return cosine; +} +} // namespace detail + /** Compute the cosine of the angle between three particles. * - * Also return all intermediate quantities: normalized vectors - * @f$ \vec{r_{ij}} @f$ (from particle @f$ j @f$ to particle @f$ i @f$) - * and @f$ \vec{r_{kj}} @f$, and their normalization constants. - * - * @param[in] box_geo Box geometry. - * @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] vec1 Vector from central particle to left particle. + * @param[in] vec2 Vector from central particle to right particle. * @param[in] sanitize_cosine Sanitize the cosine of the angle. * @return @f$ \vec{r_{ij}} @f$, @f$ \vec{r_{kj}} @f$, * @f$ \left\|\vec{r_{ij}}\right\|^{-1} @f$, * @f$ \left\|\vec{r_{kj}}\right\|^{-1} @f$, * @f$ \cos(\theta_{ijk}) @f$ */ -inline std::tuple -calc_vectors_and_cosine(BoxGeometry const &box_geo, - Utils::Vector3d const &r_mid, - Utils::Vector3d const &r_left, - Utils::Vector3d const &r_right, - bool sanitize_cosine = false) { - /* normalized vector from p_mid to p_left */ - auto vec1 = box_geo.get_mi_vector(r_left, r_mid); - auto const d1i = 1.0 / vec1.norm(); - vec1 *= d1i; - /* normalized vector from p_mid to p_right */ - auto vec2 = box_geo.get_mi_vector(r_right, r_mid); - auto const d2i = 1.0 / vec2.norm(); - vec2 *= d2i; +inline double calc_cosine(Utils::Vector3d const &vec1, + Utils::Vector3d const &vec2, + bool sanitize_cosine = false) { /* cosine of the angle between vec1 and vec2 */ - auto cosine = vec1 * vec2; + auto cos_phi = (vec1 * vec2) / std::sqrt(vec1.norm2() * vec2.norm2()); if (sanitize_cosine) { - if (cosine > TINY_COS_VALUE) - cosine = TINY_COS_VALUE; - if (cosine < -TINY_COS_VALUE) - cosine = -TINY_COS_VALUE; + cos_phi = detail::sanitize_cosine(cos_phi); } - return std::make_tuple(vec1, vec2, d1i, d2i, cosine); + return cos_phi; } /** Compute a three-body angle interaction force. @@ -77,10 +67,8 @@ calc_vectors_and_cosine(BoxGeometry const &box_geo, * See the details in @ref bondedIA_angle_force. The @f$ K(\theta_{ijk}) @f$ * term is provided as a lambda function in @p forceFactor. * - * @param[in] box_geo Box geometry. - * @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] vec1 Vector from central particle to left particle. + * @param[in] vec2 Vector from central particle to right particle. * @param[in] forceFactor Angle force term. * @param[in] sanitize_cosine Sanitize the cosine of the angle. * @tparam ForceFactor Function evaluating the angle force term @@ -89,17 +77,21 @@ calc_vectors_and_cosine(BoxGeometry const &box_geo, */ template std::tuple -angle_generic_force(BoxGeometry const &box_geo, Utils::Vector3d const &r_mid, - Utils::Vector3d const &r_left, - Utils::Vector3d const &r_right, ForceFactor forceFactor, - bool sanitize_cosine) { - auto const [vec1, vec2, d1i, d2i, cosine] = - calc_vectors_and_cosine(box_geo, r_mid, r_left, r_right, sanitize_cosine); +angle_generic_force(Utils::Vector3d const &vec1, Utils::Vector3d const &vec2, + ForceFactor forceFactor, bool sanitize_cosine) { + auto const d1 = vec1.norm(); + auto const d2 = vec2.norm(); + auto cos_phi = (vec1 * vec2) / (d1 * d2); + if (sanitize_cosine) { + cos_phi = detail::sanitize_cosine(cos_phi); + } /* force factor */ - auto const fac = forceFactor(cosine); + auto const fac = forceFactor(cos_phi); /* distribute forces */ - auto f_left = (fac * d1i) * (vec1 * cosine - vec2); - auto f_right = (fac * d2i) * (vec2 * cosine - vec1); + auto const v1 = vec1 / d1; + auto const v2 = vec2 / d2; + auto f_left = (fac / d1) * (v1 * cos_phi - v2); + auto f_right = (fac / d2) * (v2 * cos_phi - v1); auto f_mid = -(f_left + f_right); return std::make_tuple(f_mid, f_left, f_right); } diff --git a/src/core/bonded_interactions/angle_cosine.hpp b/src/core/bonded_interactions/angle_cosine.hpp index 1a75b7db7b5..4e455b0b1fc 100644 --- a/src/core/bonded_interactions/angle_cosine.hpp +++ b/src/core/bonded_interactions/angle_cosine.hpp @@ -53,11 +53,8 @@ struct AngleCosineBond { AngleCosineBond(double bend, double phi0); std::tuple - forces(BoxGeometry const &box_geo, Utils::Vector3d const &r_mid, - Utils::Vector3d const &r_left, Utils::Vector3d const &r_right) const; - double energy(BoxGeometry const &box_geo, Utils::Vector3d const &r_mid, - Utils::Vector3d const &r_left, - Utils::Vector3d const &r_right) const; + forces(Utils::Vector3d const &vec1, Utils::Vector3d const &vec2) const; + double energy(Utils::Vector3d const &vec1, Utils::Vector3d const &vec2) const; private: friend boost::serialization::access; @@ -71,17 +68,13 @@ struct AngleCosineBond { }; /** Compute the three-body angle interaction force. - * @param[in] box_geo Box geometry. - * @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] vec1 Vector from central particle to left particle. + * @param[in] vec2 Vector from central particle to right particle. * @return Forces on the second, first and third particles, in that order. */ inline std::tuple -AngleCosineBond::forces(BoxGeometry const &box_geo, - Utils::Vector3d const &r_mid, - Utils::Vector3d const &r_left, - Utils::Vector3d const &r_right) const { +AngleCosineBond::forces(Utils::Vector3d const &vec1, + Utils::Vector3d const &vec2) const { auto forceFactor = [this](double const cos_phi) { auto const sin_phi = sqrt(1 - Utils::sqr(cos_phi)); @@ -90,23 +83,16 @@ AngleCosineBond::forces(BoxGeometry const &box_geo, return -bend * (sin_phi * cos_phi0 - cos_phi * sin_phi0) / sin_phi; }; - return angle_generic_force(box_geo, r_mid, r_left, r_right, forceFactor, - false); + return angle_generic_force(vec1, vec2, forceFactor, false); } /** Computes the three-body angle interaction energy. - * @param[in] box_geo Box geometry. - * @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] vec1 Vector from central particle to left particle. + * @param[in] vec2 Vector from central particle to right particle. */ -inline double AngleCosineBond::energy(BoxGeometry const &box_geo, - Utils::Vector3d const &r_mid, - Utils::Vector3d const &r_left, - Utils::Vector3d const &r_right) const { - auto const vectors = - calc_vectors_and_cosine(box_geo, r_mid, r_left, r_right, true); - auto const cos_phi = std::get<4>(vectors); +inline double AngleCosineBond::energy(Utils::Vector3d const &vec1, + Utils::Vector3d const &vec2) const { + auto const cos_phi = calc_cosine(vec1, vec2, true); auto const sin_phi = sqrt(1 - Utils::sqr(cos_phi)); // potential: U(phi) = k * [1 - cos(phi - phi0)] // trig identity: cos(phi - phi0) = cos(phi)cos(phi0) + sin(phi)sin(phi0) diff --git a/src/core/bonded_interactions/angle_cossquare.hpp b/src/core/bonded_interactions/angle_cossquare.hpp index b21b4fd6bbf..b29d6c9cbeb 100644 --- a/src/core/bonded_interactions/angle_cossquare.hpp +++ b/src/core/bonded_interactions/angle_cossquare.hpp @@ -49,11 +49,8 @@ struct AngleCossquareBond { AngleCossquareBond(double bend, double phi0); std::tuple - forces(BoxGeometry const &box_geo, Utils::Vector3d const &r_mid, - Utils::Vector3d const &r_left, Utils::Vector3d const &r_right) const; - double energy(BoxGeometry const &box_geo, Utils::Vector3d const &r_mid, - Utils::Vector3d const &r_left, - Utils::Vector3d const &r_right) const; + forces(Utils::Vector3d const &vec1, Utils::Vector3d const &vec2) const; + double energy(Utils::Vector3d const &vec1, Utils::Vector3d const &vec2) const; private: friend boost::serialization::access; @@ -66,39 +63,28 @@ struct AngleCossquareBond { }; /** Compute the three-body angle interaction force. - * @param[in] box_geo Box geometry. - * @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] vec1 Vector from central particle to left particle. + * @param[in] vec2 Vector from central particle to right particle. * @return Forces on the second, first and third particles, in that order. */ inline std::tuple -AngleCossquareBond::forces(BoxGeometry const &box_geo, - Utils::Vector3d const &r_mid, - Utils::Vector3d const &r_left, - Utils::Vector3d const &r_right) const { +AngleCossquareBond::forces(Utils::Vector3d const &vec1, + Utils::Vector3d const &vec2) const { auto forceFactor = [this](double const cos_phi) { return bend * (cos_phi - cos_phi0); }; - return angle_generic_force(box_geo, r_mid, r_left, r_right, forceFactor, - false); + return angle_generic_force(vec1, vec2, forceFactor, false); } /** Computes the three-body angle interaction energy. - * @param[in] box_geo Box geometry. - * @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] vec1 Vector from central particle to left particle. + * @param[in] vec2 Vector from central particle to right particle. */ -inline double AngleCossquareBond::energy(BoxGeometry const &box_geo, - Utils::Vector3d const &r_mid, - Utils::Vector3d const &r_left, - Utils::Vector3d const &r_right) const { - auto const vectors = - calc_vectors_and_cosine(box_geo, r_mid, r_left, r_right, true); - auto const cos_phi = std::get<4>(vectors); +inline double AngleCossquareBond::energy(Utils::Vector3d const &vec1, + Utils::Vector3d const &vec2) const { + auto const cos_phi = calc_cosine(vec1, vec2, true); return 0.5 * bend * Utils::sqr(cos_phi - cos_phi0); } diff --git a/src/core/bonded_interactions/angle_harmonic.hpp b/src/core/bonded_interactions/angle_harmonic.hpp index e6b9e03137a..f0c813d3766 100644 --- a/src/core/bonded_interactions/angle_harmonic.hpp +++ b/src/core/bonded_interactions/angle_harmonic.hpp @@ -51,11 +51,8 @@ struct AngleHarmonicBond { } std::tuple - forces(BoxGeometry const &box_geo, Utils::Vector3d const &r_mid, - Utils::Vector3d const &r_left, Utils::Vector3d const &r_right) const; - double energy(BoxGeometry const &box_geo, Utils::Vector3d const &r_mid, - Utils::Vector3d const &r_left, - Utils::Vector3d const &r_right) const; + forces(Utils::Vector3d const &vec1, Utils::Vector3d const &vec2) const; + double energy(Utils::Vector3d const &vec1, Utils::Vector3d const &vec2) const; private: friend boost::serialization::access; @@ -67,17 +64,13 @@ struct AngleHarmonicBond { }; /** Compute the three-body angle interaction force. - * @param[in] box_geo Box geometry. - * @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] vec1 Vector from central particle to left particle. + * @param[in] vec2 Vector from central particle to right particle. * @return Forces on the second, first and third particles, in that order. */ inline std::tuple -AngleHarmonicBond::forces(BoxGeometry const &box_geo, - Utils::Vector3d const &r_mid, - Utils::Vector3d const &r_left, - Utils::Vector3d const &r_right) const { +AngleHarmonicBond::forces(Utils::Vector3d const &vec1, + Utils::Vector3d const &vec2) const { auto forceFactor = [this](double const cos_phi) { auto const sin_phi = sqrt(1 - Utils::sqr(cos_phi)); @@ -85,23 +78,16 @@ AngleHarmonicBond::forces(BoxGeometry const &box_geo, return -bend * (phi - phi0) / sin_phi; }; - return angle_generic_force(box_geo, r_mid, r_left, r_right, forceFactor, - true); + return angle_generic_force(vec1, vec2, forceFactor, true); } /** Compute the three-body angle interaction energy. - * @param[in] box_geo Box geometry. - * @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] vec1 Vector from central particle to left particle. + * @param[in] vec2 Vector from central particle to right particle. */ -inline double AngleHarmonicBond::energy(BoxGeometry const &box_geo, - Utils::Vector3d const &r_mid, - Utils::Vector3d const &r_left, - Utils::Vector3d const &r_right) const { - auto const vectors = - calc_vectors_and_cosine(box_geo, r_mid, r_left, r_right, true); - auto const cos_phi = std::get<4>(vectors); +inline double AngleHarmonicBond::energy(Utils::Vector3d const &vec1, + Utils::Vector3d const &vec2) const { + auto const cos_phi = calc_cosine(vec1, vec2, true); auto const phi = acos(cos_phi); return 0.5 * bend * Utils::sqr(phi - phi0); } diff --git a/src/core/bonded_interactions/bonded_interactions.dox b/src/core/bonded_interactions/bonded_interactions.dox index 314415a1378..b559a6c6996 100644 --- a/src/core/bonded_interactions/bonded_interactions.dox +++ b/src/core/bonded_interactions/bonded_interactions.dox @@ -134,8 +134,8 @@ * @code * # enable boost::variant with more than 20 types * target_compile_options( - * EspressoCppFlags INTERFACE -DBOOST_MPL_CFG_NO_PREPROCESSED_HEADERS - * -DBOOST_MPL_LIMIT_LIST_SIZE=40) + * espresso_cpp_flags INTERFACE -DBOOST_MPL_CFG_NO_PREPROCESSED_HEADERS + * -DBOOST_MPL_LIMIT_LIST_SIZE=40) * @endcode * * In forces_inline.hpp: * - A call to the new bond's force calculation needs to be placed in either diff --git a/src/core/bonded_interactions/bonded_tab.hpp b/src/core/bonded_interactions/bonded_tab.hpp index 515ae986340..70739fad20d 100644 --- a/src/core/bonded_interactions/bonded_tab.hpp +++ b/src/core/bonded_interactions/bonded_tab.hpp @@ -30,7 +30,6 @@ #include "config/config.hpp" -#include "BoxGeometry.hpp" #include "TabulatedPotential.hpp" #include "angle_common.hpp" #include "bonded_interactions/dihedral.hpp" @@ -93,11 +92,8 @@ struct TabulatedAngleBond : public TabulatedBond { TabulatedAngleBond(double min, double max, std::vector const &energy, std::vector const &force); std::tuple - forces(BoxGeometry const &box_geo, Utils::Vector3d const &r_mid, - Utils::Vector3d const &r_left, Utils::Vector3d const &r_right) const; - double energy(BoxGeometry const &box_geo, Utils::Vector3d const &r_mid, - Utils::Vector3d const &r_left, - Utils::Vector3d const &r_right) const; + forces(Utils::Vector3d const &vec1, Utils::Vector3d const &vec2) const; + double energy(Utils::Vector3d const &vec1, Utils::Vector3d const &vec2) const; }; /** Parameters for 4-body tabulated potential. */ @@ -111,14 +107,11 @@ struct TabulatedDihedralBond : public TabulatedBond { std::vector const &force); boost::optional> - forces(BoxGeometry const &box_geo, Utils::Vector3d const &r1, - Utils::Vector3d const &r2, Utils::Vector3d const &r3, - Utils::Vector3d const &r4) const; - boost::optional energy(BoxGeometry const &box_geo, - Utils::Vector3d const &r1, - Utils::Vector3d const &r2, - Utils::Vector3d const &r3, - Utils::Vector3d const &r4) const; + forces(Utils::Vector3d const &v12, Utils::Vector3d const &v23, + Utils::Vector3d const &v34) const; + boost::optional energy(Utils::Vector3d const &v12, + Utils::Vector3d const &v23, + Utils::Vector3d const &v34) const; }; /** Compute a tabulated bond length force. @@ -159,17 +152,13 @@ TabulatedDistanceBond::energy(Utils::Vector3d const &dx) const { } /** Compute the three-body angle interaction force. - * @param[in] box_geo Box geometry. - * @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] vec1 Vector from central particle to left particle. + * @param[in] vec2 Vector from central particle to right particle. * @return Forces on the second, first and third particles, in that order. */ inline std::tuple -TabulatedAngleBond::forces(BoxGeometry const &box_geo, - Utils::Vector3d const &r_mid, - Utils::Vector3d const &r_left, - Utils::Vector3d const &r_right) const { +TabulatedAngleBond::forces(Utils::Vector3d const &vec1, + Utils::Vector3d const &vec2) const { auto forceFactor = [this](double const cos_phi) { auto const sin_phi = sqrt(1 - Utils::sqr(cos_phi)); @@ -183,26 +172,19 @@ TabulatedAngleBond::forces(BoxGeometry const &box_geo, return -gradient / sin_phi; }; - return angle_generic_force(box_geo, r_mid, r_left, r_right, forceFactor, - true); + return angle_generic_force(vec1, vec2, forceFactor, true); } /** Compute the three-body angle interaction energy. * It is assumed that the potential is tabulated * for all angles between 0 and Pi. * - * @param[in] box_geo Box geometry. - * @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] vec1 Vector from central particle to left particle. + * @param[in] vec2 Vector from central particle to right particle. */ -inline double TabulatedAngleBond::energy(BoxGeometry const &box_geo, - Utils::Vector3d const &r_mid, - Utils::Vector3d const &r_left, - Utils::Vector3d const &r_right) const { - auto const vectors = - calc_vectors_and_cosine(box_geo, r_mid, r_left, r_right, true); - auto const cos_phi = std::get<4>(vectors); +inline double TabulatedAngleBond::energy(Utils::Vector3d const &vec1, + Utils::Vector3d const &vec2) const { + auto const cos_phi = calc_cosine(vec1, vec2, true); /* calculate phi */ #ifdef TABANGLEMINUS auto const phi = acos(-cos_phi); @@ -216,30 +198,25 @@ inline double TabulatedAngleBond::energy(BoxGeometry const &box_geo, * The forces have a singularity at @f$ \phi = 0 @f$ and @f$ \phi = \pi @f$ * (see @cite swope92a page 592). * - * @param[in] box_geo Box geometry. - * @param[in] r1 Position of the first particle. - * @param[in] r2 Position of the second particle. - * @param[in] r3 Position of the third particle. - * @param[in] r4 Position of the fourth particle. + * @param[in] v12 Vector from @p p1 to @p p2 + * @param[in] v23 Vector from @p p2 to @p p3 + * @param[in] v34 Vector from @p p3 to @p p4 * @return the forces on @p p2, @p p1, @p p3 */ inline boost::optional> -TabulatedDihedralBond::forces(BoxGeometry const &box_geo, - Utils::Vector3d const &r1, - Utils::Vector3d const &r2, - Utils::Vector3d const &r3, - Utils::Vector3d const &r4) const { +TabulatedDihedralBond::forces(Utils::Vector3d const &v12, + Utils::Vector3d const &v23, + Utils::Vector3d const &v34) const { /* vectors for dihedral angle calculation */ - Utils::Vector3d v12, v23, v34, v12Xv23, v23Xv34; + Utils::Vector3d v12Xv23, v23Xv34; double l_v12Xv23, l_v23Xv34; /* dihedral angle, cosine of the dihedral angle */ double phi, cos_phi; /* dihedral angle */ - auto const angle_is_undefined = - calc_dihedral_angle(box_geo, r1, r2, r3, r4, v12, v23, v34, v12Xv23, - l_v12Xv23, v23Xv34, l_v23Xv34, cos_phi, phi); + auto const angle_is_undefined = calc_dihedral_angle( + v12, v23, v34, v12Xv23, l_v12Xv23, v23Xv34, l_v23Xv34, cos_phi, phi); /* dihedral angle not defined - force zero */ if (angle_is_undefined) { return {}; @@ -267,24 +244,21 @@ TabulatedDihedralBond::forces(BoxGeometry const &box_geo, /** Compute the four-body dihedral interaction energy. * The energy doesn't have any singularity if the angle phi is well-defined. * - * @param[in] box_geo Box geometry. - * @param[in] r1 Position of the first particle. - * @param[in] r2 Position of the second particle. - * @param[in] r3 Position of the third particle. - * @param[in] r4 Position of the fourth particle. + * @param[in] v12 Vector from @p p1 to @p p2 + * @param[in] v23 Vector from @p p2 to @p p3 + * @param[in] v34 Vector from @p p3 to @p p4 */ -inline boost::optional TabulatedDihedralBond::energy( - BoxGeometry const &box_geo, Utils::Vector3d const &r1, - Utils::Vector3d const &r2, Utils::Vector3d const &r3, - Utils::Vector3d const &r4) const { +inline boost::optional +TabulatedDihedralBond::energy(Utils::Vector3d const &v12, + Utils::Vector3d const &v23, + Utils::Vector3d const &v34) const { /* vectors for dihedral calculations. */ - Utils::Vector3d v12, v23, v34, v12Xv23, v23Xv34; + Utils::Vector3d v12Xv23, v23Xv34; double l_v12Xv23, l_v23Xv34; /* dihedral angle, cosine of the dihedral angle */ double phi, cos_phi; - auto const angle_is_undefined = - calc_dihedral_angle(box_geo, r1, r2, r3, r4, v12, v23, v34, v12Xv23, - l_v12Xv23, v23Xv34, l_v23Xv34, cos_phi, phi); + auto const angle_is_undefined = calc_dihedral_angle( + v12, v23, v34, v12Xv23, l_v12Xv23, v23Xv34, l_v23Xv34, cos_phi, phi); /* dihedral angle not defined - energy zero */ if (angle_is_undefined) { return {}; diff --git a/src/core/bonded_interactions/dihedral.hpp b/src/core/bonded_interactions/dihedral.hpp index 13beb54c740..f5281dd0061 100644 --- a/src/core/bonded_interactions/dihedral.hpp +++ b/src/core/bonded_interactions/dihedral.hpp @@ -27,7 +27,6 @@ * the maximal bond length! */ -#include "BoxGeometry.hpp" #include "config/config.hpp" #include @@ -56,15 +55,12 @@ struct DihedralBond { boost::optional> - forces(BoxGeometry const &box_geo, Utils::Vector3d const &r1, - Utils::Vector3d const &r2, Utils::Vector3d const &r3, - Utils::Vector3d const &r4) const; + forces(Utils::Vector3d const &v12, Utils::Vector3d const &v23, + Utils::Vector3d const &v34) const; - inline boost::optional energy(BoxGeometry const &box_geo, - Utils::Vector3d const &r1, - Utils::Vector3d const &r2, - Utils::Vector3d const &r3, - Utils::Vector3d const &r4) const; + boost::optional energy(Utils::Vector3d const &v12, + Utils::Vector3d const &v23, + Utils::Vector3d const &v34) const; private: friend boost::serialization::access; @@ -86,11 +82,9 @@ struct DihedralBond { * If the a,b or b,c are parallel the dihedral angle is not defined in which * case the function returns true. Calling functions should check for that. * - * @param[in] box_geo Box geometry. - * @param[in] r1 , r2 , r3 , r4 Positions of the particles forming the dihedral - * @param[out] a Vector from @p p1 to @p p2 - * @param[out] b Vector from @p p2 to @p p3 - * @param[out] c Vector from @p p3 to @p p4 + * @param[in] a Vector from @p p1 to @p p2 + * @param[in] b Vector from @p p2 to @p p3 + * @param[in] c Vector from @p p3 to @p p4 * @param[out] aXb Vector product of a and b * @param[out] l_aXb |aXB| * @param[out] bXc Vector product of b and c @@ -99,16 +93,11 @@ struct DihedralBond { * @param[out] phi Dihedral angle in the range [0, pi] * @return Whether the angle is undefined. */ -inline bool -calc_dihedral_angle(BoxGeometry const &box_geo, Utils::Vector3d const &r1, - Utils::Vector3d const &r2, Utils::Vector3d const &r3, - Utils::Vector3d const &r4, Utils::Vector3d &a, - Utils::Vector3d &b, Utils::Vector3d &c, - Utils::Vector3d &aXb, double &l_aXb, Utils::Vector3d &bXc, - double &l_bXc, double &cosphi, double &phi) { - a = box_geo.get_mi_vector(r2, r1); - b = box_geo.get_mi_vector(r3, r2); - c = box_geo.get_mi_vector(r4, r3); +inline 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) { /* calculate vector product a X b and b X c */ aXb = vector_product(a, b); @@ -144,28 +133,24 @@ calc_dihedral_angle(BoxGeometry const &box_geo, Utils::Vector3d const &r1, * The forces have a singularity at @f$ \phi = 0 @f$ and @f$ \phi = \pi @f$ * (see @cite swope92a page 592). * - * @param[in] box_geo Box geometry. - * @param[in] r1 Position of the first particle. - * @param[in] r2 Position of the second particle. - * @param[in] r3 Position of the third particle. - * @param[in] r4 Position of the fourth particle. + * @param[in] v12 Vector from @p p1 to @p p2 + * @param[in] v23 Vector from @p p2 to @p p3 + * @param[in] v34 Vector from @p p3 to @p p4 * @return the forces on @p p2, @p p1, @p p3 */ inline boost::optional> -DihedralBond::forces(BoxGeometry const &box_geo, Utils::Vector3d const &r1, - Utils::Vector3d const &r2, Utils::Vector3d const &r3, - Utils::Vector3d const &r4) const { +DihedralBond::forces(Utils::Vector3d const &v12, Utils::Vector3d const &v23, + Utils::Vector3d const &v34) const { /* vectors for dihedral angle calculation */ - Utils::Vector3d v12, v23, v34, v12Xv23, v23Xv34; + Utils::Vector3d v12Xv23, v23Xv34; double l_v12Xv23, l_v23Xv34; /* dihedral angle, cosine of the dihedral angle */ double phi, cos_phi, sin_mphi_over_sin_phi; /* dihedral angle */ - auto const angle_is_undefined = - calc_dihedral_angle(box_geo, r1, r2, r3, r4, v12, v23, v34, v12Xv23, - l_v12Xv23, v23Xv34, l_v23Xv34, cos_phi, phi); + auto const angle_is_undefined = calc_dihedral_angle( + v12, v23, v34, v12Xv23, l_v12Xv23, v23Xv34, l_v23Xv34, cos_phi, phi); /* dihedral angle not defined - force zero */ if (angle_is_undefined) { return {}; @@ -203,25 +188,21 @@ DihedralBond::forces(BoxGeometry const &box_geo, Utils::Vector3d const &r1, /** Compute the four-body dihedral interaction energy. * The energy doesn't have any singularity if the angle phi is well-defined. * - * @param[in] box_geo Box geometry. - * @param[in] r1 Position of the first particle. - * @param[in] r2 Position of the second particle. - * @param[in] r3 Position of the third particle. - * @param[in] r4 Position of the fourth particle. + * @param[in] v12 Vector from @p p1 to @p p2 + * @param[in] v23 Vector from @p p2 to @p p3 + * @param[in] v34 Vector from @p p3 to @p p4 */ inline boost::optional -DihedralBond::energy(BoxGeometry const &box_geo, Utils::Vector3d const &r1, - Utils::Vector3d const &r2, Utils::Vector3d const &r3, - Utils::Vector3d const &r4) const { +DihedralBond::energy(Utils::Vector3d const &v12, Utils::Vector3d const &v23, + Utils::Vector3d const &v34) const { /* vectors for dihedral calculations. */ - Utils::Vector3d v12, v23, v34, v12Xv23, v23Xv34; + Utils::Vector3d v12Xv23, v23Xv34; double l_v12Xv23, l_v23Xv34; /* dihedral angle, cosine of the dihedral angle */ double phi, cos_phi; - auto const angle_is_undefined = - calc_dihedral_angle(box_geo, r1, r2, r3, r4, v12, v23, v34, v12Xv23, - l_v12Xv23, v23Xv34, l_v23Xv34, cos_phi, phi); + auto const angle_is_undefined = calc_dihedral_angle( + v12, v23, v34, v12Xv23, l_v12Xv23, v23Xv34, l_v23Xv34, cos_phi, phi); /* dihedral angle not defined - energy zero */ if (angle_is_undefined) { return {}; diff --git a/src/core/energy_inline.hpp b/src/core/energy_inline.hpp index f70bf8d578e..c693051f83d 100644 --- a/src/core/energy_inline.hpp +++ b/src/core/energy_inline.hpp @@ -246,17 +246,19 @@ calc_bonded_energy(Bonded_IA_Parameters const &iaparams, Particle const &p1, throw BondUnknownTypeError(); } // 1 partner if (n_partners == 2) { + auto const vec1 = box_geo.get_mi_vector(p2->pos(), p1.pos()); + auto const vec2 = box_geo.get_mi_vector(p3->pos(), p1.pos()); if (auto const *iap = boost::get(&iaparams)) { - return iap->energy(box_geo, p1.pos(), p2->pos(), p3->pos()); + return iap->energy(vec1, vec2); } if (auto const *iap = boost::get(&iaparams)) { - return iap->energy(box_geo, p1.pos(), p2->pos(), p3->pos()); + return iap->energy(vec1, vec2); } if (auto const *iap = boost::get(&iaparams)) { - return iap->energy(box_geo, p1.pos(), p2->pos(), p3->pos()); + return iap->energy(vec1, vec2); } if (auto const *iap = boost::get(&iaparams)) { - return iap->energy(box_geo, p1.pos(), p2->pos(), p3->pos()); + return iap->energy(vec1, vec2); } if (boost::get(&iaparams)) { runtimeWarningMsg() << "Unsupported bond type " + @@ -267,11 +269,15 @@ calc_bonded_energy(Bonded_IA_Parameters const &iaparams, Particle const &p1, throw BondUnknownTypeError(); } // 2 partners if (n_partners == 3) { + // note: particles in a dihedral bond are ordered as p2-p1-p3-p4 + auto const v12 = box_geo.get_mi_vector(p1.pos(), p2->pos()); + auto const v23 = box_geo.get_mi_vector(p3->pos(), p1.pos()); + auto const v34 = box_geo.get_mi_vector(p4->pos(), p3->pos()); if (auto const *iap = boost::get(&iaparams)) { - return iap->energy(box_geo, p2->pos(), p1.pos(), p3->pos(), p4->pos()); + return iap->energy(v12, v23, v34); } if (auto const *iap = boost::get(&iaparams)) { - return iap->energy(box_geo, p2->pos(), p1.pos(), p3->pos(), p4->pos()); + return iap->energy(v12, v23, v34); } if (boost::get(&iaparams)) { runtimeWarningMsg() << "Unsupported bond type " + diff --git a/src/core/forces_inline.hpp b/src/core/forces_inline.hpp index 387ce2fe3d3..2d21f2483be 100644 --- a/src/core/forces_inline.hpp +++ b/src/core/forces_inline.hpp @@ -359,22 +359,24 @@ inline boost::optional< calc_bonded_three_body_force(Bonded_IA_Parameters const &iaparams, BoxGeometry const &box_geo, Particle const &p1, Particle const &p2, Particle const &p3) { + auto const vec1 = box_geo.get_mi_vector(p2.pos(), p1.pos()); + auto const vec2 = box_geo.get_mi_vector(p3.pos(), p1.pos()); if (auto const *iap = boost::get(&iaparams)) { - return iap->forces(box_geo, p1.pos(), p2.pos(), p3.pos()); + return iap->forces(vec1, vec2); } if (auto const *iap = boost::get(&iaparams)) { - return iap->forces(box_geo, p1.pos(), p2.pos(), p3.pos()); + return iap->forces(vec1, vec2); } if (auto const *iap = boost::get(&iaparams)) { - return iap->forces(box_geo, p1.pos(), p2.pos(), p3.pos()); + return iap->forces(vec1, vec2); } #ifdef TABULATED if (auto const *iap = boost::get(&iaparams)) { - return iap->forces(box_geo, p1.pos(), p2.pos(), p3.pos()); + return iap->forces(vec1, vec2); } #endif if (auto const *iap = boost::get(&iaparams)) { - return iap->calc_forces(box_geo, p1, p2, p3); + return iap->calc_forces(vec1, vec2); } throw BondUnknownTypeError(); } @@ -412,12 +414,16 @@ calc_bonded_four_body_force(Bonded_IA_Parameters const &iaparams, if (auto const *iap = boost::get(&iaparams)) { return iap->calc_forces(box_geo, p1, p2, p3, p4); } + // note: particles in a dihedral bond are ordered as p2-p1-p3-p4 + auto const v12 = box_geo.get_mi_vector(p1.pos(), p2.pos()); + auto const v23 = box_geo.get_mi_vector(p3.pos(), p1.pos()); + auto const v34 = box_geo.get_mi_vector(p4.pos(), p3.pos()); if (auto const *iap = boost::get(&iaparams)) { - return iap->forces(box_geo, p2.pos(), p1.pos(), p3.pos(), p4.pos()); + return iap->forces(v12, v23, v34); } #ifdef TABULATED if (auto const *iap = boost::get(&iaparams)) { - return iap->forces(box_geo, p2.pos(), p1.pos(), p3.pos(), p4.pos()); + return iap->forces(v12, v23, v34); } #endif throw BondUnknownTypeError(); diff --git a/src/core/immersed_boundary/ibm_triel.cpp b/src/core/immersed_boundary/ibm_triel.cpp index 19d512924cd..b5c80cb4c07 100644 --- a/src/core/immersed_boundary/ibm_triel.cpp +++ b/src/core/immersed_boundary/ibm_triel.cpp @@ -20,7 +20,6 @@ #include "immersed_boundary/ibm_triel.hpp" #include "BoxGeometry.hpp" -#include "Particle.hpp" #include "ibm_common.hpp" #include "system/System.hpp" @@ -73,17 +72,15 @@ RotateForces(Utils::Vector2d const &f1_rot, Utils::Vector2d const &f2_rot, } // namespace boost::optional> -IBMTriel::calc_forces(BoxGeometry const &box_geo, Particle const &p1, - Particle const &p2, Particle const &p3) const { +IBMTriel::calc_forces(Utils::Vector3d const &vec1, + Utils::Vector3d const &vec2) const { // Calculate the current shape of the triangle (l,lp,cos(phi),sin(phi)); // l = length between 1 and 3 // get_mi_vector is an ESPResSo function which considers PBC - auto const vec2 = box_geo.get_mi_vector(p3.pos(), p1.pos()); auto const l = vec2.norm(); // lp = length between 1 and 2 - auto const vec1 = box_geo.get_mi_vector(p2.pos(), p1.pos()); auto const lp = vec1.norm(); // Check for sanity @@ -93,7 +90,7 @@ IBMTriel::calc_forces(BoxGeometry const &box_geo, Particle const &p1, // angles between these vectors; calculated directly via the products auto const cosPhi = (vec1 * vec2) / (lp * l); - auto const vecpro = vector_product(vec1, vec2); + auto const vecpro = Utils::vector_product(vec1, vec2); auto const sinPhi = vecpro.norm() / (l * lp); // Displacement gradient tensor D: eq. (C.9) in @cite kruger12a @@ -217,7 +214,7 @@ IBMTriel::IBMTriel(const int ind1, const int ind2, const int ind3, // cospo / sinpo angle functions between these vectors; calculated directly // via the products cosPhi0 = (temp_l0 * temp_lp0) / (l0 * lp0); - auto const vecpro = vector_product(temp_l0, temp_lp0); + auto const vecpro = Utils::vector_product(temp_l0, temp_lp0); sinPhi0 = vecpro.norm() / (l0 * lp0); // Use the values determined above for further constants of the stretch-force diff --git a/src/core/immersed_boundary/ibm_triel.hpp b/src/core/immersed_boundary/ibm_triel.hpp index 78c52027a67..429ce876a44 100644 --- a/src/core/immersed_boundary/ibm_triel.hpp +++ b/src/core/immersed_boundary/ibm_triel.hpp @@ -20,8 +20,6 @@ #ifndef IBM_TRIEL_H #define IBM_TRIEL_H -#include "BoxGeometry.hpp" -#include "Particle.hpp" #include "config/config.hpp" #include @@ -70,8 +68,7 @@ struct IBMTriel { * @return the forces on @p p1, @p p2, @p p3 */ boost::optional> - calc_forces(BoxGeometry const &box_geo, Particle const &p1, - Particle const &p2, Particle const &p3) const; + calc_forces(Utils::Vector3d const &vec1, Utils::Vector3d const &vec2) const; private: friend boost::serialization::access;