/*! \brief harmonic bond type
*
* It represents the interaction of the form
- * V(r, forceConstant, equilDistance) = 0.5 * forceConstant * (r - equilConstant)^2
+ * V(r, forceConstant, equilConstant) = 0.5 * forceConstant * (r - equilConstant)^2
*/
using HarmonicBondType = TwoParameterInteraction<struct HarmonicBondTypeParameter>;
/*! \brief GROMOS bond type
*
* It represents the interaction of the form
- * V(r, forceConstant, equilDistance) = 0.25 * forceConstant * (r^2 - equilConstant^2)^2
+ * V(r, forceConstant, equilConstant) = 0.25 * forceConstant * (r^2 - equilConstant^2)^2
*/
class G96BondType : public TwoParameterInteraction<struct G96BondTypeParameter>
{
/*! \brief FENE bond type
*
* It represents the interaction of the form
- * V(r, forceConstant, equilDistance) = - 0.5 * forceConstant * equilDistance^2 * log( 1 - (r / equilConstant)^2)
+ * V(r, forceConstant, equilConstant) = - 0.5 * forceConstant * equilConstant^2 * log( 1 - (r / equilConstant)^2)
*/
using FENEBondType = TwoParameterInteraction<struct FENEBondTypeParameter>;
/*! \brief Half-attractive quartic bond type
*
* It represents the interaction of the form
- * V(r, forceConstant, equilConstant) = 0.5 * forceConstant * (r - equilDistance)^4
+ * V(r, forceConstant, equilConstant) = 0.5 * forceConstant * (r - equilConstant)^4
*/
using HalfAttractiveQuarticBondType =
TwoParameterInteraction<struct HalfAttractiveQuarticBondTypeParameter>;
*
* It represents the interaction of the form
* V(r, quadraticForceConstant, cubicForceConstant, equilConstant) = quadraticForceConstant * (r -
- * equilDistance)^2 + quadraticForceConstant * cubicForceConstant * (r - equilConstant)
+ * equilConstant)^2 + quadraticForceConstant * cubicForceConstant * (r - equilConstant)
*/
class CubicBondType
{
/*! \brief Morse bond type
*
* It represents the interaction of the form
- * V(r, forceConstant, exponent, equilDistance) = forceConstant * ( 1 - exp( -exponent * (r - equilConstant))
+ * V(r, forceConstant, exponent, equilConstant) = forceConstant * ( 1 - exp( -exponent * (r - equilConstant))
*/
class MorseBondType
{
/*! \brief Harmonic angle type
*
* It represents the interaction of the form
- * V(theta, forceConstant, equilAngle) = 0.5 * forceConstant * (theta - equilAngle)^2
+ * V(theta, forceConstant, equilConstant) = 0.5 * forceConstant * (theta - equilConstant)^2
*/
using HarmonicAngle = AngleInteractionType<struct HarmonicAngleParameter>;
/*! \brief G96 or Cosine-based angle type
*
* This represents the interaction of the form
- * V(cos(theta), forceConstant, cos(equilAngle)) = 0.5 * forceConstant * (cos(theta) - cos(equilAngle))^2
+ * V(cos(theta), forceConstant, cos(equilConstant)) = 0.5 * forceConstant * (cos(theta) - cos(equilConstant))^2
*/
using G96Angle = CosineParamAngle<struct G96AngleParameter>;
/*! \brief Restricted angle type
*
* This represents the interaction of the form
- * V(cos(theta), forceConstant, cos(equilAngle)) =
- * 0.5 * forceConstant * (cos(theta) - cos(equilAngle))^2 / (sin(theta))^2
+ * V(cos(theta), forceConstant, cos(equilConstant)) =
+ * 0.5 * forceConstant * (cos(theta) - cos(equilConstant))^2 / (sin(theta))^2
*/
using RestrictedAngle = CosineParamAngle<struct RestrictedAngleParameter>;
/*! \brief Quartic angle type
*
* It represents the interaction of the form of a fourth order polynomial
- * V(theta, forceConstant, equilAngle) = sum[i = 0 -> 4](forceConstant_i * (theta - equilAngle)^i
+ * V(theta, forceConstant, equilConstant) = sum[i = 0 -> 4](forceConstant_i * (theta - equilConstant)^i
*/
class QuarticAngle
{
}
//! \brief dummy op in case shift forces are not computed (will be removed by the compiler)
-inline std::nullptr_t accessShiftForces([[maybe_unused]] int shiftIndex,
- [[maybe_unused]] gmx::ArrayRef<std::nullptr_t> shiftForces)
+inline std::nullptr_t accessShiftForces(int /* shiftIndex */,
+ gmx::ArrayRef<std::nullptr_t> /* shiftForces */)
{
return nullptr;
}
BasicVector* fi, BasicVector* fj, BasicVector* fk,
ShiftForce* shf_ij, ShiftForce* shf_kj, ShiftForce* shf_c)
{
- if (parameters.manifest == ThreeCenterType::Cargo::ij)
+if (parameters.manifest == ThreeCenterType::Cargo::ij)
{
// i-j bond
return computeTwoCenter(parameters.twoCenter(), rij, fi, fj, shf_ij, shf_c);
template<class ThreeCenterType, class BasicVector, class ShiftForce>
inline NBLIB_ALWAYS_INLINE
std::enable_if_t<!HasTwoCenterAggregate<ThreeCenterType>::value, BasicVectorValueType_t<BasicVector>>
-addTwoCenterAggregate([[maybe_unused]] const ThreeCenterType& parameters,
- [[maybe_unused]] const BasicVector& rij,
- [[maybe_unused]] const BasicVector& rkj,
- [[maybe_unused]] BasicVector* fi,
- [[maybe_unused]] BasicVector* fj,
- [[maybe_unused]] BasicVector* fk,
- [[maybe_unused]] ShiftForce* shf_ij,
- [[maybe_unused]] ShiftForce* shf_kj,
- [[maybe_unused]] ShiftForce* shf_c)
+addTwoCenterAggregate(const ThreeCenterType& /* parameters */,
+ const BasicVector& /* rij */,
+ const BasicVector& /* rkj */,
+ BasicVector* /* fi */,
+ BasicVector* /* fj */,
+ BasicVector* /* fk */,
+ ShiftForce* /* shf_ij */,
+ ShiftForce* /* shf_kj */,
+ ShiftForce* /* shf_c */)
{
return 0.0;
};
{
using ValueType = BasicVectorValueType_t<BasicVector>;
// 42 flops from the norm() calls
- auto [ci, cj, ck, energy] = crossBondAngleKernel(norm(rij), norm(rkj), norm(rik), parameters);
+ auto [ci, cj, ck, energy] = threeCenterKernel(norm(rij), norm(rkj), norm(rik), parameters);
BasicVector fi_loc = ci * rij + ck * rik;
BasicVector fk_loc = cj * rkj - ck * rik;
template<class FourCenterType, class BasicVector>
inline NBLIB_ALWAYS_INLINE
std::enable_if_t<!HasThreeCenterAggregate<FourCenterType>::value, BasicVectorValueType_t<BasicVector>>
-addThreeCenterAggregate([[maybe_unused]] const FourCenterType& parameters,
- [[maybe_unused]] const BasicVector& rij,
- [[maybe_unused]] const BasicVector& rkj,
- [[maybe_unused]] const BasicVector& rkl,
- [[maybe_unused]] BasicVector* fi,
- [[maybe_unused]] BasicVector* fj,
- [[maybe_unused]] BasicVector* fk,
- [[maybe_unused]] BasicVector* fl)
+addThreeCenterAggregate(const FourCenterType& /* parameters*/,
+ const BasicVector& /* rij */,
+ const BasicVector& /* rkj */,
+ const BasicVector& /* rkl */,
+ BasicVector* /* fi */,
+ BasicVector* /* fj */,
+ BasicVector* /* fk */,
+ BasicVector* /* fl */)
{
return 0.0;
};
FiveCenterType fiveCenterTypeParams = parameters[std::get<5>(index)];
// this dispatch function is not in use yet, because CMap is not yet implemented
- // we don't want to add [[maybe_unused]] in the signature
- // and we also don't want compiler warnings, so we cast to void
+ // we don't want to add [[maybe_unused]] in the signature since the params will
+ // be used once CMap is implemented, and we also don't want compiler warnings,
+ // so we cast to void.
(void)fiveCenterTypeParams;
(void)forces;
#include "gromacs/math/vectypes.h"
#include "nblib/listed_forces/bondtypes.h"
-#define NBLIB_ALWAYS_INLINE __attribute((always_inline))
-
namespace nblib
{
/*! \brief kernel to calculate the scalar part for linear angle forces
* for lambda = 0
*
- * \param k force constant
- * \param a contribution of rij vector
- * \param dr weighted distance metric
+ * \param k force constant
+ * \param a0 equilibrium angle
+ * \param angle current angle vaule
*
* \return tuple<force, potential energy>
*/
template <class T>
-inline std::tuple<T, T, T> linearAnglesScalar(T k, T a, T dr)
+inline std::tuple<T, T, T> linearAnglesScalarForce(T k, T a0, T angle)
{
- T b = T(1.0) - a;
+ T b = T(1.0) - a0;
- T kdr = k * dr;
- T epot = 0.5 * kdr * dr;
+ T kdr = k * angle;
+ T epot = 0.5 * kdr * angle;
- T ci = a * k;
+ T ci = a0 * k;
T ck = b * k;
return std::make_tuple(ci, ck, epot);
template <class T>
inline auto threeCenterKernel(T dr, const LinearAngle& angle)
{
- return linearAnglesScalar(angle.forceConstant(), angle.equilConstant(), dr);
+ return linearAnglesScalarForce(angle.forceConstant(), angle.equilConstant(), dr);
}
//! Harmonic Angle
* \param rkj input bond length between particles k & j
* \param rik input bond length between particles i & k
*
- * \return tuple<force, potential energy>
+ * \return tuple<atom i force, atom j force, atom k force, potential energy>
*/
template <class T>
T epot = k * sik * (sij + skj);
- T ci = -k * sik / rij;
- T cj = -k * sik / rkj;
- T ck = -k * (sij + skj) / rik;
+ T fi = -k * sik / rij;
+ T fj = -k * sik / rkj;
+ T fk = -k * (sij + skj) / rik;
- return std::make_tuple(ci, cj, ck, epot);
+ return std::make_tuple(fi, fj, fk, epot);
/* That was 13 flops */
}
//! Cross bond-bond interaction
template <class T>
-inline auto crossBondAngleKernel(T drij, T drkj, T drik, const CrossBondAngle& crossBondAngle)
+inline auto threeCenterKernel(T drij, T drkj, T drik, const CrossBondAngle& crossBondAngle)
{
return crossBondAngleScalarForce(crossBondAngle.forceConstant(), crossBondAngle.equilDistanceIJ(), crossBondAngle.equilDistanceKJ(), crossBondAngle.equilDistanceIK(), drij, drkj, drik);
}
{
}
-/*! \brief Spreads and accumulates the bonded forces to the two atoms and adds the virial contribution when needed
+/*! \brief Spreads and accumulates the forces between two atoms and adds the virial contribution when needed
*
- * \p shiftIndex is used as the periodic shift.
+ * \tparam T The type of vector, e.g. float, double, etc
+ * \param force The Force
+ * \param dx Distance between centers
+ * \param force_i Force on center i
+ * \param force_j Force on center j
+ * \param shf_ik Shift force between centers i and j
+ * \param shf_c Shift force at the "center" of the two center interaction
*/
template <class T, class ShiftForce>
-inline void spreadTwoCenterForces(const T bondForce,
+inline void spreadTwoCenterForces(const T force,
const gmx::BasicVector<T>& dx,
gmx::BasicVector<T>* force_i,
gmx::BasicVector<T>* force_j,
ShiftForce* shf_ij,
ShiftForce* shf_c)
{
- gmx::BasicVector<T> fij = bondForce * dx;
+ gmx::BasicVector<T> fij = force * dx;
*force_i += fij;
*force_j -= fij;
/*! \brief spread force to 3 centers based on scalar force and angle
*
- * \tparam T
- * \param cos_theta
- * \param force
- * \param r_ij
- * \param r_kj
- * \param force_i
- * \param force_j
- * \param force_k
+ * \tparam T The type of vector, e.g. float, double, etc
+ * \param cos_theta Angle between two vectors
+ * \param force The Force
+ * \param dxIJ Distance between centers i and j
+ * \param dxJK Distance between centers j and k
+ * \param force_i Force on center i
+ * \param force_j Force on center j
+ * \param force_k Force on center k
+ * \param shf_ik Shift force between centers i and j
+ * \param shf_kj Shift force between centers k and j
+ * \param shf_c Shift force at the center subtending the angle
*/
template <class T, class ShiftForce>
inline void spreadThreeCenterForces(T cos_theta,
T force,
- const gmx::BasicVector<T>& r_ij,
- const gmx::BasicVector<T>& r_kj,
+ const gmx::BasicVector<T>& dxIJ,
+ const gmx::BasicVector<T>& dxKJ,
gmx::BasicVector<T>* force_i,
gmx::BasicVector<T>* force_j,
gmx::BasicVector<T>* force_k,
{
T st = force / std::sqrt(1 - cos_theta2); /* 12 */
T sth = st * cos_theta; /* 1 */
- T nrij2 = dot(r_ij, r_ij); /* 5 */
- T nrkj2 = dot(r_kj, r_kj); /* 5 */
+ T nrij2 = dot(dxIJ, dxIJ); /* 5 */
+ T nrkj2 = dot(dxKJ, dxKJ); /* 5 */
T cik = st / std::sqrt(nrij2 * nrkj2); /* 11 */
T cii = sth / nrij2; /* 1 */
T ckk = sth / nrkj2; /* 1 */
/* 33 */
- gmx::BasicVector<T> f_i = cii * r_ij - cik * r_kj;
- gmx::BasicVector<T> f_k = ckk * r_kj - cik * r_ij;
+ gmx::BasicVector<T> f_i = cii * dxIJ - cik * dxKJ;
+ gmx::BasicVector<T> f_k = ckk * dxKJ - cik * dxIJ;
gmx::BasicVector<T> f_j = T(-1.0) * (f_i + f_k);
*force_i += f_i;
*force_j += f_j;
}
//! Four-center category common
+
+/*! \brief spread force to 4 centers
+ *
+ * \tparam T The type of vector, e.g. float, double, etc
+ * \param dxIJ Distance between centers i and j
+ * \param dxKJ Distance between centers j and k
+ * \param dxKL Distance between centers k and l
+ * \param m Cross product of \p dxIJ x \p dxKJ
+ * \param m Cross product of \p dxKJ x \p dxKL
+ * \param force_i Force on center i
+ * \param force_j Force on center j
+ * \param force_k Force on center k
+ * \param force_k Force on center l
+ * \param shf_ik Shift force between centers i and j
+ * \param shf_kj Shift force between centers k and j
+ * \param shf_lj Shift force between centers k and j
+ * \param shf_c Shift force at the center subtending the angle
+ */
template <class T, class ShiftForce>
inline void spreadFourCenterForces(T force,
const gmx::BasicVector<T>& dxIJ,
} // namespace nblib
-#undef NBLIB_ALWAYS_INLINE
-
#endif // NBLIB_LISTEDFORCES_KERNELS_HPP
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff