Unverified Commit 3387ef0d authored by bors[bot]'s avatar bors[bot] Committed by GitHub

[ci skip][skip ci][skip netlify] -bors-staging-tmp-3100

parents 3f5c608d 2a32c8e3
Pipeline #8924 skipped
......@@ -529,10 +529,6 @@ section :ref:`Isotropic non-bonded interactions`):
Some of the short-range interactions have additional features:
- ``LJ_WARN_WHEN_CLOSE`` This adds an additional check to the Lennard-Jones
potentials that prints a warning if particles come so close to each other
that the simulation becomes unphysical.
- ``LJGEN_SOFTCORE`` This modifies the generic Lennard-Jones potential
(``LENNARD_JONES_GENERIC``) with tunable parameters.
......
......@@ -64,9 +64,8 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
#define VIRTUAL_SITES_RELATIVE
#define EXPERIMENTAL_FEATURES
// DEBUG Switches
#define LJ_WARN_WHEN_CLOSE
// DEBUG Switches
#define COMM_DEBUG
#define EVENT_DEBUG
#define CELL_DEBUG
......
......@@ -71,7 +71,6 @@ EK_DOUBLE_PREC requires ELECTROKINETICS
TABULATED
LENNARD_JONES
WCA
LJ_WARN_WHEN_CLOSE
LENNARD_JONES_GENERIC implies LENNARD_JONES
LJCOS
LJCOS2
......
......@@ -29,7 +29,6 @@
#include "angle_common.hpp"
#include "bonded_interaction_data.hpp"
#include "grid.hpp"
#include "particle_data.hpp"
#include <utils/math/sqr.hpp>
......@@ -39,16 +38,16 @@
int angle_cosine_set_params(int bond_type, double bend, double phi0);
/** Compute the three-body angle interaction force.
* @param[in] p_mid Second/middle particle.
* @param[in] p_left First/left particle.
* @param[in] p_right Third/right particle.
* @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] iaparams Bonded parameters for the angle interaction.
* @return Forces on the second, first and third particles, in that order.
*/
inline std::tuple<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>
calc_angle_cosine_3body_forces(Particle const *const p_mid,
Particle const *const p_left,
Particle const *const p_right,
calc_angle_cosine_3body_forces(Utils::Vector3d const &r_mid,
Utils::Vector3d const &r_left,
Utils::Vector3d const &r_right,
Bonded_ia_parameters const *const iaparams) {
auto forceFactor = [&iaparams](double const cos_phi) {
......@@ -61,14 +60,13 @@ calc_angle_cosine_3body_forces(Particle const *const p_mid,
return -k * (sin_phi * cos_phi0 - cos_phi * sin_phi0) / sin_phi;
};
return calc_angle_generic_force(p_mid->r.p, p_left->r.p, p_right->r.p,
forceFactor, false);
return calc_angle_generic_force(r_mid, r_left, r_right, forceFactor, false);
}
/** Compute the three-body angle interaction force.
* @param[in] p_mid Second/middle particle.
* @param[in] p_left First/left particle.
* @param[in] p_right Third/right particle.
* @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] iaparams Bonded parameters for the angle interaction.
* @param[out] f_mid Force on @p p_mid.
* @param[out] f_left Force on @p p_left.
......@@ -76,29 +74,28 @@ calc_angle_cosine_3body_forces(Particle const *const p_mid,
* @retval false
*/
inline bool calc_angle_cosine_force(
Particle const *const p_mid, Particle const *const p_left,
Particle const *const p_right, Bonded_ia_parameters const *const iaparams,
Utils::Vector3d const &r_mid, Utils::Vector3d const &r_left,
Utils::Vector3d const &r_right, Bonded_ia_parameters const *const iaparams,
Utils::Vector3d &f_mid, Utils::Vector3d &f_left, Utils::Vector3d &f_right) {
std::tie(f_mid, f_left, f_right) =
calc_angle_cosine_3body_forces(p_mid, p_left, p_right, iaparams);
calc_angle_cosine_3body_forces(r_mid, r_left, r_right, iaparams);
return false;
}
/** Computes the three-body angle interaction energy.
* @param[in] p_mid Second/middle particle.
* @param[in] p_left First/left particle.
* @param[in] p_right Third/right particle.
* @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] iaparams Bonded parameters for the angle interaction.
* @param[out] _energy Energy.
* @retval false
*/
inline bool angle_cosine_energy(Particle const *const p_mid,
Particle const *const p_left,
Particle const *const p_right,
inline bool angle_cosine_energy(Utils::Vector3d const &r_mid,
Utils::Vector3d const &r_left,
Utils::Vector3d const &r_right,
Bonded_ia_parameters const *const iaparams,
double *_energy) {
auto const vectors =
calc_vectors_and_cosine(p_mid->r.p, p_left->r.p, p_right->r.p, true);
auto const vectors = calc_vectors_and_cosine(r_mid, r_left, r_right, true);
auto const cos_phi = std::get<4>(vectors);
auto const sin_phi = sqrt(1 - Utils::sqr(cos_phi));
auto const cos_phi0 = iaparams->p.angle_cosine.cos_phi0;
......
......@@ -29,7 +29,6 @@
#include "angle_common.hpp"
#include "bonded_interaction_data.hpp"
#include "grid.hpp"
#include "particle_data.hpp"
#include <utils/constants.hpp>
#include <utils/math/sqr.hpp>
......@@ -40,16 +39,16 @@
int angle_cossquare_set_params(int bond_type, double bend, double phi0);
/** Compute the three-body angle interaction force.
* @param[in] p_mid Second/middle particle.
* @param[in] p_left First/left particle.
* @param[in] p_right Third/right particle.
* @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] iaparams Bonded parameters for the angle interaction.
* @return Forces on the second, first and third particles, in that order.
*/
inline std::tuple<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>
calc_angle_cossquare_3body_forces(Particle const *const p_mid,
Particle const *const p_left,
Particle const *const p_right,
calc_angle_cossquare_3body_forces(Utils::Vector3d const &r_mid,
Utils::Vector3d const &r_left,
Utils::Vector3d const &r_right,
Bonded_ia_parameters const *const iaparams) {
auto forceFactor = [&iaparams](double const cos_phi) {
......@@ -58,14 +57,13 @@ calc_angle_cossquare_3body_forces(Particle const *const p_mid,
return k * (cos_phi - cos_phi0);
};
return calc_angle_generic_force(p_mid->r.p, p_left->r.p, p_right->r.p,
forceFactor, false);
return calc_angle_generic_force(r_mid, r_left, r_right, forceFactor, false);
}
/** Compute the three-body angle interaction force.
* @param[in] p_mid Second/middle particle.
* @param[in] p_left First/left particle.
* @param[in] p_right Third/right particle.
* @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] iaparams Bonded parameters for the angle interaction.
* @param[out] f_mid Force on @p p_mid.
* @param[out] f_left Force on @p p_left.
......@@ -73,29 +71,28 @@ calc_angle_cossquare_3body_forces(Particle const *const p_mid,
* @retval false
*/
inline bool calc_angle_cossquare_force(
Particle const *const p_mid, Particle const *const p_left,
Particle const *const p_right, Bonded_ia_parameters const *const iaparams,
Utils::Vector3d const &r_mid, Utils::Vector3d const &r_left,
Utils::Vector3d const &r_right, Bonded_ia_parameters const *const iaparams,
Utils::Vector3d &f_mid, Utils::Vector3d &f_left, Utils::Vector3d &f_right) {
std::tie(f_mid, f_left, f_right) =
calc_angle_cossquare_3body_forces(p_mid, p_left, p_right, iaparams);
calc_angle_cossquare_3body_forces(r_mid, r_left, r_right, iaparams);
return false;
}
/** Computes the three-body angle interaction energy.
* @param[in] p_mid Second/middle particle.
* @param[in] p_left First/left particle.
* @param[in] p_right Third/right particle.
* @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] iaparams Bonded parameters for the angle interaction.
* @param[out] _energy Energy.
* @retval false
*/
inline bool angle_cossquare_energy(Particle const *const p_mid,
Particle const *const p_left,
Particle const *const p_right,
inline bool angle_cossquare_energy(Utils::Vector3d const &r_mid,
Utils::Vector3d const &r_left,
Utils::Vector3d const &r_right,
Bonded_ia_parameters const *const iaparams,
double *_energy) {
auto const vectors =
calc_vectors_and_cosine(p_mid->r.p, p_left->r.p, p_right->r.p, true);
auto const vectors = calc_vectors_and_cosine(r_mid, r_left, r_right, true);
auto const cos_phi = std::get<4>(vectors);
auto const cos_phi0 = iaparams->p.angle_cossquare.cos_phi0;
auto const k = iaparams->p.angle_cossquare.bend;
......
......@@ -26,10 +26,8 @@
* @ref bondedIA_angle_harmonic.
*/
#include "bonded_interaction_data.hpp"
#include "particle_data.hpp"
#include "angle_common.hpp"
#include "bonded_interaction_data.hpp"
#include "grid.hpp"
#include <utils/math/sqr.hpp>
......@@ -40,16 +38,16 @@
int angle_harmonic_set_params(int bond_type, double bend, double phi0);
/** Compute the three-body angle interaction force.
* @param p_mid Second/middle particle.
* @param p_left First/left particle.
* @param p_right Third/right particle.
* @param iaparams Bonded parameters for the angle interaction.
* @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] iaparams Bonded parameters for the angle interaction.
* @return Forces on the second, first and third particles, in that order.
*/
inline std::tuple<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>
calc_angle_harmonic_3body_forces(Particle const *const p_mid,
Particle const *const p_left,
Particle const *const p_right,
calc_angle_harmonic_3body_forces(Utils::Vector3d const &r_mid,
Utils::Vector3d const &r_left,
Utils::Vector3d const &r_right,
Bonded_ia_parameters const *const iaparams) {
auto forceFactor = [&iaparams](double const cos_phi) {
......@@ -60,14 +58,13 @@ calc_angle_harmonic_3body_forces(Particle const *const p_mid,
return -k * (phi - phi0) / sin_phi;
};
return calc_angle_generic_force(p_mid->r.p, p_left->r.p, p_right->r.p,
forceFactor, true);
return calc_angle_generic_force(r_mid, r_left, r_right, forceFactor, true);
}
/** Compute the three-body angle interaction force.
* @param[in] p_mid Second/middle particle.
* @param[in] p_left First/left particle.
* @param[in] p_right Third/right particle.
* @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] iaparams Bonded parameters for the angle interaction.
* @param[out] f_mid Force on @p p_mid.
* @param[out] f_left Force on @p p_left.
......@@ -75,29 +72,28 @@ calc_angle_harmonic_3body_forces(Particle const *const p_mid,
* @retval false
*/
inline bool calc_angle_harmonic_force(
Particle const *const p_mid, Particle const *const p_left,
Particle const *const p_right, Bonded_ia_parameters const *const iaparams,
Utils::Vector3d const &r_mid, Utils::Vector3d const &r_left,
Utils::Vector3d const &r_right, Bonded_ia_parameters const *const iaparams,
Utils::Vector3d &f_mid, Utils::Vector3d &f_left, Utils::Vector3d &f_right) {
std::tie(f_mid, f_left, f_right) =
calc_angle_harmonic_3body_forces(p_mid, p_left, p_right, iaparams);
calc_angle_harmonic_3body_forces(r_mid, r_left, r_right, iaparams);
return false;
}
/** Compute the three-body angle interaction energy.
* @param[in] p_mid Second/middle particle.
* @param[in] p_left First/left particle.
* @param[in] p_right Third/right particle.
* @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] iaparams Bonded parameters for the angle interaction.
* @param[out] _energy Energy.
* @retval false
*/
inline bool angle_harmonic_energy(Particle const *const p_mid,
Particle const *const p_left,
Particle const *const p_right,
inline bool angle_harmonic_energy(Utils::Vector3d const &r_mid,
Utils::Vector3d const &r_left,
Utils::Vector3d const &r_right,
Bonded_ia_parameters const *const iaparams,
double *_energy) {
auto const vectors =
calc_vectors_and_cosine(p_mid->r.p, p_left->r.p, p_right->r.p, true);
auto const vectors = calc_vectors_and_cosine(r_mid, r_left, r_right, true);
auto const cos_phi = std::get<4>(vectors);
auto const phi = acos(cos_phi);
auto const phi0 = iaparams->p.angle_harmonic.phi0;
......
......@@ -27,15 +27,11 @@
* Implementation in \ref bonded_coulomb.cpp
*/
/************************************************************/
#include "config.hpp"
#ifdef ELECTROSTATICS
#include "bonded_interaction_data.hpp"
#include "debug.hpp"
#include "particle_data.hpp"
/** Set the parameters for the bonded Coulomb potential
*
......@@ -45,42 +41,36 @@
int bonded_coulomb_set_params(int bond_type, double prefactor);
/** Compute the bonded Coulomb pair force.
* @param[in] p1 First particle.
* @param[in] p2 Second particle.
* @param[in] q1q2 Product of the particle charges.
* @param[in] iaparams Interaction parameters.
* @param[in] dx %Distance between the particles.
* @param[out] force Force.
* @retval false
*/
inline bool calc_bonded_coulomb_pair_force(
Particle const *const p1, Particle const *const p2,
Bonded_ia_parameters const *const iaparams, Utils::Vector3d const &dx,
Utils::Vector3d &force) {
double const q1q2, Bonded_ia_parameters const *const iaparams,
Utils::Vector3d const &dx, Utils::Vector3d &force) {
auto const dist2 = dx.norm2();
auto const dist = std::sqrt(dist2);
auto const fac =
iaparams->p.bonded_coulomb.prefactor * p1->p.q * p2->p.q / (dist * dist2);
auto const dist3 = dist2 * std::sqrt(dist2);
auto const fac = iaparams->p.bonded_coulomb.prefactor * q1q2 / dist3;
force = fac * dx;
return false;
}
/** Compute the bonded Coulomb pair energy.
* @param[in] p1 First particle.
* @param[in] p2 Second particle.
* @param[in] q1q2 Product of the particle charges.
* @param[in] iaparams Interaction parameters.
* @param[in] dx %Distance between the particles.
* @param[out] _energy Energy.
* @retval false
*/
inline bool
bonded_coulomb_pair_energy(Particle const *const p1, Particle const *const p2,
bonded_coulomb_pair_energy(double const q1q2,
Bonded_ia_parameters const *const iaparams,
Utils::Vector3d const &dx, double *_energy) {
auto const dist = dx.norm();
*_energy = iaparams->p.bonded_coulomb.prefactor * p1->p.q * p2->p.q / dist;
*_energy = iaparams->p.bonded_coulomb.prefactor * q1q2 / dist;
return false;
}
......
......@@ -28,14 +28,11 @@
* Implementation in \ref bonded_coulomb_sr.cpp.
*/
/************************************************************/
#include "config.hpp"
#ifdef ELECTROSTATICS
#include "bonded_interaction_data.hpp"
#include "debug.hpp"
#include "electrostatics_magnetostatics/coulomb_inline.hpp"
#include "particle_data.hpp"
......
......@@ -61,14 +61,10 @@
*
* * The recommended signatures of the force and energy functions are:
* @code{.cpp}
* inline bool calc_fene_pair_force(Particle const *const p1,
* Particle const *const p2,
* Bonded_ia_parameters const *const iaparams,
* inline bool calc_fene_pair_force(Bonded_ia_parameters const *const iaparams,
* Vector3d const &dx, Vector3d &force);
* inline bool fene_pair_energy(Particle const *const p1,
* Particle const *const p2,
* Bonded_ia_parameters const *const iaparams,
* Vector3d const &dx, double *energy);
* inline bool fene_pair_energy(Bonded_ia_parameters const *const iaparams,
* Vector3d const &dx, double *_energy);
* @endcode
* * The setter function gets a \p bond_type which is a numerical id
* identifying the number of the bond type in the simulation.
......@@ -125,7 +121,7 @@
* code looks like this:
* @code{.cpp}
* case BONDED_IA_FENE:
* bond_broken = calc_fene_pair_force(p1, p2, iaparams, dx, force);
* bond_broken = calc_fene_pair_force(iaparams, dx, force);
* break;
* @endcode
* - in function @ref calc_bond_pair_force(): add the new interaction to the
......@@ -133,8 +129,8 @@
* the code looks like this:
* @code{.cpp}
* case BONDED_IA_ANGLE_HARMONIC:
* bond_broken = calc_angle_harmonic_force(p1, p2, p3, iaparams, force,
* force2, force3);
* bond_broken = calc_angle_harmonic_force(p1->r.p, p2->r.p, p3->r.p,
* iaparams, force, force2, force3);
* break;
* @endcode
* * In energy_inline.hpp:
......@@ -143,7 +139,7 @@
* switch statement. For the FENE bond, e.g., the code looks like this:
* @code{.cpp}
* case BONDED_IA_FENE:
* bond_broken = fene_pair_energy(p1, p2, iaparams, dx, &ret);
* bond_broken = fene_pair_energy(iaparams, dx, &ret);
* break;
* @endcode
* * Pressure, stress tensor and virial calculation: if your bonded
......
......@@ -33,8 +33,6 @@
#include "angle_common.hpp"
#include "bonded_interactions/bonded_interaction_data.hpp"
#include "bonded_interactions/dihedral.hpp"
#include "debug.hpp"
#include "particle_data.hpp"
#include <tuple>
#include <utils/constants.hpp>
......@@ -107,16 +105,16 @@ inline bool tab_bond_energy(Bonded_ia_parameters const *const iaparams,
}
/** Compute the three-body angle interaction force.
* @param p_mid Second/middle particle.
* @param p_left First/left particle.
* @param p_right Third/right particle.
* @param iaparams Bonded parameters for the angle interaction.
* @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] iaparams Bonded parameters for the angle interaction.
* @return Forces on the second, first and third particles, in that order.
*/
inline std::tuple<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>
calc_angle_3body_tabulated_forces(Particle const *const p_mid,
Particle const *const p_left,
Particle const *const p_right,
calc_angle_3body_tabulated_forces(Utils::Vector3d const &r_mid,
Utils::Vector3d const &r_left,
Utils::Vector3d const &r_right,
Bonded_ia_parameters const *const iaparams) {
auto forceFactor = [&iaparams](double const cos_phi) {
......@@ -131,14 +129,13 @@ calc_angle_3body_tabulated_forces(Particle const *const p_mid,
return -gradient / sin_phi;
};
return calc_angle_generic_force(p_mid->r.p, p_left->r.p, p_right->r.p,
forceFactor, true);
return calc_angle_generic_force(r_mid, r_left, r_right, forceFactor, true);
}
/** Compute the three-body angle interaction force.
* @param[in] p_mid Second/middle particle.
* @param[in] p_left First/left particle.
* @param[in] p_right Third/right particle.
* @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] iaparams Bonded parameters for the angle interaction.
* @param[out] f_mid Force on @p p_mid.
* @param[out] f_left Force on @p p_left.
......@@ -146,12 +143,12 @@ calc_angle_3body_tabulated_forces(Particle const *const p_mid,
* @retval false
*/
inline bool calc_tab_angle_force(
Particle const *const p_mid, Particle const *const p_left,
Particle const *const p_right, Bonded_ia_parameters const *const iaparams,
Utils::Vector3d const &r_mid, Utils::Vector3d const &r_left,
Utils::Vector3d const &r_right, Bonded_ia_parameters const *const iaparams,
Utils::Vector3d &f_mid, Utils::Vector3d &f_left, Utils::Vector3d &f_right) {
std::tie(f_mid, f_left, f_right) =
calc_angle_3body_tabulated_forces(p_mid, p_left, p_right, iaparams);
calc_angle_3body_tabulated_forces(r_mid, r_left, r_right, iaparams);
return false;
}
......@@ -159,20 +156,19 @@ inline bool calc_tab_angle_force(
* It is assumed that the potential is tabulated
* for all angles between 0 and Pi.
*
* @param[in] p_mid Second/middle particle.
* @param[in] p_left First/left particle.
* @param[in] p_right Third/right particle.
* @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] iaparams Bonded parameters for the angle interaction.
* @param[out] _energy Energy.
* @retval false
*/
inline bool tab_angle_energy(Particle const *const p_mid,
Particle const *const p_left,
Particle const *const p_right,
inline bool tab_angle_energy(Utils::Vector3d const &r_mid,
Utils::Vector3d const &r_left,
Utils::Vector3d const &r_right,
Bonded_ia_parameters const *const iaparams,
double *_energy) {
auto const vectors =
calc_vectors_and_cosine(p_mid->r.p, p_left->r.p, p_right->r.p, true);
auto const vectors = calc_vectors_and_cosine(r_mid, r_left, r_right, true);
auto const cos_phi = std::get<4>(vectors);
/* calculate phi */
#ifdef TABANGLEMINUS
......@@ -187,10 +183,10 @@ inline bool tab_angle_energy(Particle const *const p_mid,
/** Compute the four-body dihedral interaction force.
* This function is not tested yet.
*
* @param[in] p2 Second particle.
* @param[in] p1 First particle.
* @param[in] p3 Third particle.
* @param[in] p4 Fourth particle.
* @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] iaparams Bonded parameters for the dihedral interaction.
* @param[out] force2 Force on particle 2.
* @param[out] force1 Force on particle 1.
......@@ -198,8 +194,8 @@ inline bool tab_angle_energy(Particle const *const p_mid,
* @return false
*/
inline bool
calc_tab_dihedral_force(Particle const *const p2, Particle const *const p1,
Particle const *const p3, Particle const *const p4,
calc_tab_dihedral_force(Utils::Vector3d const &r1, Utils::Vector3d const &r2,
Utils::Vector3d const &r3, Utils::Vector3d const &r4,
Bonded_ia_parameters const *const iaparams,
Utils::Vector3d &force2, Utils::Vector3d &force1,
Utils::Vector3d &force3) {
......@@ -212,7 +208,7 @@ calc_tab_dihedral_force(Particle const *const p2, Particle const *const p1,
auto const *tab_pot = iaparams->p.tab.pot;
/* dihedral angle */
calc_dihedral_angle(p1, p2, p3, p4, v12, v23, v34, v12Xv23, &l_v12Xv23,
calc_dihedral_angle(r1, r2, r3, r4, v12, v23, v34, v12Xv23, &l_v12Xv23,
v23Xv34, &l_v23Xv34, &cos_phi, &phi);
/* dihedral angle not defined - force zero */
if (phi == -1.0) {
......@@ -245,18 +241,18 @@ calc_tab_dihedral_force(Particle const *const p2, Particle const *const p1,
/** Compute the four-body dihedral interaction energy.
* This function is not tested yet.
*
* @param[in] p2 Second particle.
* @param[in] p1 First particle.
* @param[in] p3 Third particle.
* @param[in] p4 Fourth particle.
* @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] iaparams Bonded parameters for the dihedral interaction.
* @param[out] _energy Energy.
* @return false
*/
inline bool tab_dihedral_energy(Particle const *const p2,
Particle const *const p1,
Particle const *const p3,
Particle const *const p4,
inline bool tab_dihedral_energy(Utils::Vector3d const &r1,
Utils::Vector3d const &r2,
Utils::Vector3d const &r3,
Utils::Vector3d const &r4,
Bonded_ia_parameters const *const iaparams,
double *_energy) {
/* vectors for dihedral calculations. */
......@@ -265,7 +261,7 @@ inline bool tab_dihedral_energy(Particle const *const p2,
/* dihedral angle, cosine of the dihedral angle */
double phi, cos_phi;
auto const *tab_pot = iaparams->p.tab.pot;
calc_dihedral_angle(p1, p2, p3, p4, v12, v23, v34, v12Xv23, &l_v12Xv23,
calc_dihedral_angle(r1, r2, r3, r4, v12, v23, v34, v12Xv23, &l_v12Xv23,
v23Xv34, &l_v23Xv34, &cos_phi, &phi);
*_energy = tab_pot->energy(phi);
......
......@@ -52,7 +52,7 @@ int dihedral_set_params(int bond_type, int mult, double bend, double phase);
* If the a,b or b,c are parallel the dihedral angle is not defined in which
* case the routine returns phi=-1. Calling functions should check for that
*
* @param[in] p1 , p2 , p3 , p4 Particles forming the dihedral
* @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
......@@ -64,14 +64,14 @@ int dihedral_set_params(int bond_type, int mult, double bend, double phase);
* @param[out] phi Dihedral angle
*/
inline void
calc_dihedral_angle(Particle const *const p1, Particle const *const p2,
Particle const *const p3, Particle const *const p4,
calc_dihedral_angle(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 = get_mi_vector(p2->r.p, p1->r.p, box_geo);
b = get_mi_vector(p3->r.p, p2->r.p, box_geo);
c = get_mi_vector(p4->r.p, p3->r.p, box_geo);
a = get_mi_vector(r2, r1, box_geo);
b = get_mi_vector(r3, r2, box_geo);
c = get_mi_vector(r4, r3, box_geo);
/* calculate vector product a X b and b X c */
aXb = vector_product(a, b);
......@@ -104,10 +104,10 @@ calc_dihedral_angle(Particle const *const p1, Particle const *const p2,
/** Compute the four-body dihedral interaction force.
*
* @param[in] p2 Second particle.
* @param[in] p1 First particle.
* @param[in] p3 Third particle.
* @param[in] p4 Fourth particle.
* @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] iaparams Bonded parameters for the dihedral interaction.
* @param[out] force2 Force on particle 2.
* @param[out] force1 Force on particle 1.
......@@ -115,8 +115,8 @@ calc_dihedral_angle(Particle const *const p1, Particle const *const p2,
* @return false
*/
inline bool
calc_dihedral_force(Particle const *const p2, Particle const *const p1,
Particle const *const p3, Particle const *const p4,
calc_dihedral_force(Utils::Vector3d const &r1, Utils::Vector3d const &r2,
Utils::Vector3d const &r3, Utils::Vector3d const &r4,
Bonded_ia_parameters const *const iaparams,
Utils::Vector3d &force2, Utils::Vector3d &force1,
Utils::Vector3d &force3) {
......@@ -129,7 +129,7 @@ calc_dihedral_force(Particle const *const p2, Particle const *const p1,
double fac;
/* dihedral angle */
calc_dihedral_angle(p1, p2, p3, p4, v12, v23, v34, v12Xv23, &l_v12Xv23,
calc_dihedral_angle(r1, r2, r3, r4, v12, v23, v34, v12Xv23, &l_v12Xv23,
v23Xv34, &l_v23Xv34, &cosphi, &phi);
/* dihedral angle not defined - force zero */
if (phi == -1.0) {
......@@ -177,25 +177,25 @@ calc_dihedral_force(Particle const *const p2, Particle const *const p1,
/** Compute the four-body dihedral interaction energy.
*
* @param[in] p2 Second particle.
* @param[in] p1 First particle.
* @param[in] p3 Third particle.
* @param[in] p4 Fourth particle.
* @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] iaparams Bonded parameters for the dihedral interaction.
* @param[out] _energy Energy.
* @return false
*/
inline bool dihedral_energy(Particle const *const p1, Particle const *const p2,
Particle const *const p3, Particle const *const p4,
Bonded_ia_parameters const *const iaparams,
double *_energy) {
inline bool
dihedral_energy(Utils::Vector3d const &r1, Utils::Vector3d const &r2,