Skip to content

Commit

Permalink
core: Move BoxGeometry out of bonded IA kernels
Browse files Browse the repository at this point in the history
  • Loading branch information
@jngrad
jngrad committed Sep 20, 2023
1 parent 11e1df1 commit e945b9a
Show file tree
Hide file tree
Showing 11 changed files with 167 additions and 256 deletions.
76 changes: 34 additions & 42 deletions src/core/bonded_interactions/angle_common.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,63 +24,51 @@
* Common code for functions calculating angle forces.
*/

#include "BoxGeometry.hpp"
#include "config/config.hpp"

#include <utils/Vector.hpp>

#include <cmath>
#include <tuple>

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<Utils::Vector3d, Utils::Vector3d, double, double, double>
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.
*
* 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
Expand All @@ -89,17 +77,21 @@ calc_vectors_and_cosine(BoxGeometry const &box_geo,
*/
template <typename ForceFactor>
std::tuple<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>
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);
}
Expand Down
38 changes: 12 additions & 26 deletions src/core/bonded_interactions/angle_cosine.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,11 +53,8 @@ struct AngleCosineBond {
AngleCosineBond(double bend, double phi0);

std::tuple<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>
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;
Expand All @@ -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<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>
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));
Expand All @@ -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)
Expand Down
38 changes: 12 additions & 26 deletions src/core/bonded_interactions/angle_cossquare.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,11 +49,8 @@ struct AngleCossquareBond {
AngleCossquareBond(double bend, double phi0);

std::tuple<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>
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;
Expand All @@ -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<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>
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);
}

Expand Down
38 changes: 12 additions & 26 deletions src/core/bonded_interactions/angle_harmonic.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,11 +51,8 @@ struct AngleHarmonicBond {
}

std::tuple<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>
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;
Expand All @@ -67,41 +64,30 @@ 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<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>
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));
auto const phi = acos(cos_phi);
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);
}
Expand Down
4 changes: 2 additions & 2 deletions src/core/bonded_interactions/bonded_interactions.dox
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading

0 comments on commit e945b9a

Please sign in to comment.