namespace nblib
{
-template<class TwoCenterType, class BasicVector>
+//! \brief returns the address of an element in the shiftForces buffer
+template<class ShiftForce>
+inline ShiftForce* accessShiftForces(int shiftIndex, gmx::ArrayRef<ShiftForce> shiftForces)
+{
+ return shiftForces.data() + shiftIndex;
+}
+
+//! \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)
+{
+ return nullptr;
+}
+
+template<class TwoCenterType, class BasicVector, class ShiftForce>
inline NBLIB_ALWAYS_INLINE
-auto computeTwoCenter(const TwoCenterType& parameters, const BasicVector& dx, BasicVector* fi, BasicVector* fj)
+auto computeTwoCenter(const TwoCenterType& parameters, const BasicVector& dx, BasicVector* fi, BasicVector* fj,
+ ShiftForce* sh_f, ShiftForce* sh_fc)
{
using ValueType = BasicVectorValueType_t<BasicVector>;
if (dr2 != 0.0)
{
force /= dr;
- spreadTwoCenterForces(force, dx, fi, fj);
+ spreadTwoCenterForces(force, dx, fi, fj, sh_f, sh_fc);
}
return energy;
* \param[in] pbc Object used for computing distances accounting for PBC's
* \return Computed kernel energies
*/
-template <class Buffer, class TwoCenterType, class BasicVector, class Pbc,
+template <class Buffer, class TwoCenterType, class BasicVector, class ShiftForce, class Pbc,
std::enable_if_t<Contains<TwoCenterType, SupportedTwoCenterTypes>{}>* = nullptr>
inline NBLIB_ALWAYS_INLINE
auto dispatchInteraction(InteractionIndex<TwoCenterType> index,
gmx::ArrayRef<const TwoCenterType> bondInstances,
gmx::ArrayRef<const BasicVector> x,
Buffer* forces,
+ gmx::ArrayRef<ShiftForce> shiftForces,
const Pbc& pbc)
{
KernelEnergy<BasicVectorValueType_t<BasicVector>> energy;
int i = std::get<0>(index);
int j = std::get<1>(index);
- BasicVector x1 = x[i];
- BasicVector x2 = x[j];
+ BasicVector xi = x[i];
+ BasicVector xj = x[j];
TwoCenterType bond = bondInstances[std::get<2>(index)];
BasicVector dx;
- // calculate x1 - x2 modulo pbc
- pbc.dxAiuc(x1, x2, dx);
+ // calculate xi - xj modulo pbc
+ int sIdx = pbc.dxAiuc(xi, xj, dx);
- energy.carrier() = computeTwoCenter(bond, dx, &(*forces)[i], &(*forces)[j]);
+ ShiftForce* sh_f = accessShiftForces(sIdx, shiftForces);
+ ShiftForce* sh_fc = accessShiftForces(gmx::c_centralShiftIndex, shiftForces);
+
+ energy.carrier() = computeTwoCenter(bond, dx, &(*forces)[i], &(*forces)[j], sh_f, sh_fc);
return energy;
}
-template<class ThreeCenterType, class BasicVector>
+template<class ThreeCenterType, class BasicVector, class ShiftForce>
inline NBLIB_ALWAYS_INLINE
std::enable_if_t<HasTwoCenterAggregate<ThreeCenterType>::value, BasicVectorValueType_t<BasicVector>>
addTwoCenterAggregate(const ThreeCenterType& parameters, const BasicVector& rij, const BasicVector& rkj,
- BasicVector* fi, BasicVector* fj, BasicVector* fk)
+ BasicVector* fi, BasicVector* fj, BasicVector* fk,
+ ShiftForce* shf_ij, ShiftForce* shf_kj, ShiftForce* shf_c)
{
if (parameters.manifest == ThreeCenterType::Cargo::ij)
{
// i-j bond
- return computeTwoCenter(parameters.twoCenter(), rij, fi, fj);
+ return computeTwoCenter(parameters.twoCenter(), rij, fi, fj, shf_ij, shf_c);
}
if (parameters.manifest == ThreeCenterType::Cargo::jk)
{
// j-k bond
- return computeTwoCenter(parameters.twoCenter(), rkj, fk, fj);
+ return computeTwoCenter(parameters.twoCenter(), rkj, fk, fj, shf_kj, shf_c);
}
// aggregate is empty
return 0.0;
};
-template<class ThreeCenterType, class BasicVector>
+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& rkj,
[[maybe_unused]] BasicVector* fi,
[[maybe_unused]] BasicVector* fj,
- [[maybe_unused]] BasicVector* fk)
+ [[maybe_unused]] BasicVector* fk,
+ [[maybe_unused]] ShiftForce* shf_ij,
+ [[maybe_unused]] ShiftForce* shf_kj,
+ [[maybe_unused]] ShiftForce* shf_c)
{
return 0.0;
};
-template<class ThreeCenterType, class BasicVector>
+template<class ThreeCenterType, class BasicVector, class ShiftForce>
inline NBLIB_ALWAYS_INLINE
auto computeThreeCenter(const ThreeCenterType& parameters, const BasicVector& rij, const BasicVector& rkj,
- [[maybe_unused]] const BasicVector& rik, BasicVector* fi, BasicVector* fj, BasicVector* fk)
+ const BasicVector& /* rik */, BasicVector* fi, BasicVector* fj, BasicVector* fk,
+ ShiftForce* shf_ij, ShiftForce* shf_kj, ShiftForce* shf_c)
{
using ValueType = BasicVectorValueType_t<BasicVector>;
// calculate 3-center common quantities: angle between x1-x2 and x2-x3
ValueType costh = basicVectorCosAngle(rij, rkj); /* 25 */
ValueType theta = std::acos(costh); /* 10 */
- // call type-specific angle kernel, e.g. harmonic, linear, quartic, etc.
+ // call type-specific angle kernel, e.g. harmonic, restricted, quartic, etc.
auto [force, energy] = threeCenterKernel(theta, parameters);
- spreadThreeCenterForces(costh, force, rij, rkj, fi, fj, fk);
+ spreadThreeCenterForces(costh, force, rij, rkj, fi, fj, fk, shf_ij, shf_kj, shf_c);
return energy;
}
* \param[in] PBC
* \return Computed kernel energies
*/
-template <class Buffer, class ThreeCenterType, class BasicVector, class Pbc,
+template <class Buffer, class ThreeCenterType, class BasicVector, class ShiftForce, class Pbc,
std::enable_if_t<Contains<ThreeCenterType, SupportedThreeCenterTypes>{}>* = nullptr>
inline NBLIB_ALWAYS_INLINE
auto dispatchInteraction(InteractionIndex<ThreeCenterType> index,
gmx::ArrayRef<const ThreeCenterType> parameters,
gmx::ArrayRef<const BasicVector> x,
Buffer* forces,
+ gmx::ArrayRef<ShiftForce> shiftForces,
const Pbc& pbc)
{
KernelEnergy<BasicVectorValueType_t<BasicVector>> energy;
BasicVector fi{0,0,0}, fj{0,0,0}, fk{0,0,0};
BasicVector rij, rkj, rik;
- pbc.dxAiuc(xi, xj, rij); /* 3 */
- pbc.dxAiuc(xk, xj, rkj); /* 3 */
+ int sIdx_ij = pbc.dxAiuc(xi, xj, rij); /* 3 */
+ int sIdx_kj = pbc.dxAiuc(xk, xj, rkj); /* 3 */
pbc.dxAiuc(xi, xk, rik); /* 3 */
- energy.carrier() = computeThreeCenter(threeCenterParameters, rij, rkj, rik, &fi, &fj, &fk);
- energy.twoCenterAggregate() = addTwoCenterAggregate(threeCenterParameters, rij, rkj, &fi, &fj, &fk);
+ ShiftForce* shf_ij = accessShiftForces(sIdx_ij, shiftForces);
+ ShiftForce* shf_kj = accessShiftForces(sIdx_kj, shiftForces);
+ ShiftForce* shf_c = accessShiftForces(gmx::c_centralShiftIndex, shiftForces);
+
+ energy.carrier() = computeThreeCenter(threeCenterParameters, rij, rkj, rik, &fi, &fj, &fk, shf_ij, shf_kj, shf_c);
+ energy.twoCenterAggregate() = addTwoCenterAggregate(threeCenterParameters, rij, rkj, &fi, &fj, &fk, shf_ij, shf_kj, shf_c);
(*forces)[i] += fi;
(*forces)[j] += fj;
* \param[in] pbc Object used for computing distances accounting for PBC's
* \return Computed kernel energies
*/
-template <class Buffer, class FourCenterType, class BasicVector, class Pbc,
+template <class Buffer, class FourCenterType, class BasicVector, class ShiftForce, class Pbc,
std::enable_if_t<Contains<FourCenterType, SupportedFourCenterTypes>{}>* = nullptr>
inline NBLIB_ALWAYS_INLINE
auto dispatchInteraction(InteractionIndex<FourCenterType> index,
gmx::ArrayRef<const FourCenterType> parameters,
gmx::ArrayRef<const BasicVector> x,
Buffer* forces,
+ gmx::ArrayRef<ShiftForce> shiftForces,
const Pbc& pbc)
{
using RealScalar = BasicVectorValueType_t<BasicVector>;
BasicVector fi{0,0,0}, fj{0,0,0}, fk{0,0,0}, fl{0,0,0};
- BasicVector dxIJ, dxKJ, dxKL;
- pbc.dxAiuc(xi, xj, dxIJ);
- pbc.dxAiuc(xk, xj, dxKJ);
+ BasicVector dxIJ, dxKJ, dxKL, dxLJ;
+ int sIdx_ij = pbc.dxAiuc(xi, xj, dxIJ);
+ int sIdx_kj = pbc.dxAiuc(xk, xj, dxKJ);
pbc.dxAiuc(xk, xl, dxKL);
+ int sIdx_lj = pbc.dxAiuc(xl, xj, dxLJ);
+
+ ShiftForce* shf_ij = accessShiftForces(sIdx_ij, shiftForces);
+ ShiftForce* shf_kj = accessShiftForces(sIdx_kj, shiftForces);
+ ShiftForce* shf_lj = accessShiftForces(sIdx_lj, shiftForces);
+ ShiftForce* shf_c = accessShiftForces(gmx::c_centralShiftIndex, shiftForces);
+
FourCenterType fourCenterTypeParams = parameters[std::get<4>(index)];
BasicVector m, n;
- RealScalar phi = dihedralPhi(dxIJ, dxKJ, dxKL, m, n);
+ RealScalar phi = dihedralPhi(dxIJ, dxKJ, dxKL, &m, &n);
auto [force, kernelEnergy] = fourCenterKernel(phi, fourCenterTypeParams);
energy.carrier() = kernelEnergy;
energy.threeCenterAggregate() = addThreeCenterAggregate(fourCenterTypeParams, dxIJ, dxKJ, dxKL, &fi, &fj, &fk, &fl);
- spreadFourCenterForces(force, dxIJ, dxKJ, dxKL, m, n, &fi, &fj, &fk, &fl);
+ spreadFourCenterForces(force, dxIJ, dxKJ, dxKL, m, n, &fi, &fj, &fk, &fl, shf_ij, shf_kj, shf_lj, shf_c);
(*forces)[i] += fi;
(*forces)[j] += fj;
* \param[in] pbc Object used for computing distances accounting for PBC's
* \return Computed kernel energies
*/
-template <class Buffer, class FiveCenterType, class BasicVector, class Pbc,
+template <class Buffer, class FiveCenterType, class BasicVector, class ShiftForce, class Pbc,
std::enable_if_t<Contains<FiveCenterType, SupportedFiveCenterTypes>{}>* = nullptr>
inline NBLIB_ALWAYS_INLINE
auto dispatchInteraction(InteractionIndex<FiveCenterType> index,
gmx::ArrayRef<const FiveCenterType> parameters,
gmx::ArrayRef<const BasicVector> x,
Buffer* forces,
+ [[maybe_unused]] gmx::ArrayRef<ShiftForce> shiftForces,
const Pbc& pbc)
{
KernelEnergy<BasicVectorValueType_t<BasicVector>> energy;
* \param[in] pbc Object used for computing distances accounting for PBC's
* \return Computed kernel energies
*/
-template <class Index, class InteractionType, class Buffer, class Pbc>
+template <class Index, class InteractionType, class Buffer, class ShiftForce, class Pbc>
auto computeForces(gmx::ArrayRef<const Index> indices,
gmx::ArrayRef<const InteractionType> parameters,
gmx::ArrayRef<const Vec3> x,
Buffer* forces,
+ gmx::ArrayRef<ShiftForce> shiftForces,
const Pbc& pbc)
{
KernelEnergy<BasicVectorValueType_t<Vec3>> energy;
for (const auto& index : indices)
{
- energy += dispatchInteraction(index, parameters, x, forces, pbc);
+ energy += dispatchInteraction(index, parameters, x, forces, shiftForces, pbc);
}
return energy;
}
+//! \brief convenience overload without shift forces
+template <class Index, class InteractionType, class Buffer, class Pbc>
+auto computeForces(gmx::ArrayRef<const Index> indices,
+ gmx::ArrayRef<const InteractionType> parameters,
+ gmx::ArrayRef<const Vec3> x,
+ Buffer* forces,
+ const Pbc& pbc)
+{
+ return computeForces(indices, parameters, x, forces, gmx::ArrayRef<std::nullptr_t>{}, pbc);
+}
+
/*! \brief implement a loop over bond types and accumulate their force contributions
*
* \param[in] interactions interaction pairs and bond parameters
* \param[in] pbc Object used for computing distances accounting for PBC's
* \return Computed kernel energies
*/
-template<class Buffer, class Pbc>
+template<class Buffer, class ShiftForce, class Pbc>
auto reduceListedForces(const ListedInteractionData& interactions,
gmx::ArrayRef<const Vec3> x,
Buffer* forces,
+ gmx::ArrayRef<ShiftForce> shiftForces,
const Pbc& pbc)
{
using ValueType = BasicVectorValueType_t<Vec3>;
energies.fill(0);
// calculate one bond type
- auto computeForceType = [forces, &x, &energies, &pbc](const auto& interactionElement) {
+ auto computeForceType = [forces, x, shiftForces, &energies, &pbc](const auto& interactionElement) {
using InteractionType = typename std::decay_t<decltype(interactionElement)>::type;
gmx::ArrayRef<const InteractionIndex<InteractionType>> indices(interactionElement.indices);
gmx::ArrayRef<const InteractionType> parameters(interactionElement.parameters);
- KernelEnergy<ValueType> energy = computeForces(indices, parameters, x, forces, pbc);
+
+ KernelEnergy<ValueType> energy = computeForces(indices, parameters, x, forces, shiftForces, pbc);
energies[CarrierIndex<InteractionType>{}] += energy.carrier();
energies[TwoCenterAggregateIndex<InteractionType>{}] += energy.twoCenterAggregate();
return energies;
}
+//! \brief convenience overload without shift forces
+template<class Buffer, class Pbc>
+auto reduceListedForces(const ListedInteractionData& interactions,
+ gmx::ArrayRef<const Vec3> x,
+ Buffer* forces,
+ const Pbc& pbc)
+{
+ return reduceListedForces(interactions, x, forces, gmx::ArrayRef<std::nullptr_t>{}, pbc);
+}
+
} // namespace nblib
#undef NBLIB_ALWAYS_INLINE
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2020,2021, by the GROMACS development team, led by
+ * Copyright (c) 2020, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include <tuple>
-#include "gromacs/math/vec.h"
-#include "nblib/basicdefinitions.h"
-#include "bondtypes.h"
+#include "gromacs/math/functions.h"
+#include "gromacs/math/vectypes.h"
+#include "nblib/listed_forces/bondtypes.h"
+
+#define NBLIB_ALWAYS_INLINE __attribute((always_inline))
namespace nblib
{
T dx = x - x0;
T dx2 = dx * dx;
- T force = -k * dx * x;
- T epot = 0.25 * k * dx2;
+ T force = -k * dx;
+ T epot = 0.5 * k * dx2;
return std::make_tuple(force, epot);
- /* That was 7 flops */
+ /* That was 6 flops */
}
/*! \brief kernel to calculate the scalar part for forth power pontential bond forces
T dx = x - x0;
T dx2 = dx * dx;
- T force = -kk * dx * x;
- T epot = 0.25 * kk * dx2;
+ T force = -kk * dx;
+ T epot = 0.5 * kk * dx2;
// TODO: Check if this is 1/2 or 1/4
T dvdlambda = 0.5 * (kB - kA) * dx2 + (xA - xB) * kk * dx;
/* That was 21 flops */
}
-//! Abstraction layer for different 2-center bonds. Forth power case
+//! Abstraction layer for different 2-center bonds. Fourth power case
template <class T>
inline auto bondKernel(T dr, const G96BondType& bond)
{
- // NOTE: Not assuming GROMACS' convention of storing squared bond.equilConstant() for this type
- return g96ScalarForce(bond.forceConstant(), bond.equilConstant() * bond.equilConstant(), dr * dr);
+ auto [force, ePot] = g96ScalarForce(bond.forceConstant(), bond.equilConstant(), dr*dr);
+ force *= dr;
+ ePot *= 0.5;
+ return std::make_tuple(force, ePot);
}
T kexp = k * omexp; /* 1 */
T epot = kexp * omexp; /* 1 */
- T force = -2.0 * beta * exponent * omexp; /* 4 */
+ T force = -2.0 * beta * exponent * kexp; /* 4 */
return std::make_tuple(force, epot);
- /* That was 23 flops */
+ /* That was 20 flops */
}
/*! \brief kernel to calculate the scalar part for morse potential bond forces
T kexp = k * omexp; /* 1 */
T epot = kexp * omexp; /* 1 */
- T force = -2.0 * beta * exponent * omexp; /* 4 */
+ T force = -2.0 * beta * exponent * kexp; /* 4 */
T dvdlambda = (kB - kA) * omexp * omexp
- (2.0 - 2.0 * omexp) * omexp * k
T omx2_ox02 = 1.0 - (x2 / x02);
T epot = -0.5 * k * x02 * std::log(omx2_ox02);
- T force = -k * x / omx2_ox02;
+ T force = -k / omx2_ox02;
return std::make_tuple(force, epot);
template <class T>
inline auto bondKernel(T dr, const FENEBondType& bond)
{
- return FENEScalarForce(bond.forceConstant(), bond.equilConstant(), dr);
+ auto [force, ePot] = FENEScalarForce(bond.forceConstant(), bond.equilConstant(), dr);
+ force *= dr;
+ return std::make_tuple(force, ePot);
}
inline std::tuple<T, T> cubicScalarForce(T kc, T kq, T x0, T x)
{
T dx = x - x0;
+ //T dx2 = dx * dx;
T kdist = kq * dx;
T kdist2 = kdist * dx;
//! Three-center interaction type kernels
-// linear angles go here
-// quartic angles go here
-
-//! Three-center interaction type dispatch
+//! Harmonic Angle
template <class T>
inline auto threeCenterKernel(T dr, const HarmonicAngle& angle)
{
template <class T>
inline auto fourCenterKernel(T phi, const ProperDihedral& properDihedral)
{
- const T deltaPhi = properDihedral.multiplicity() * phi - properDihedral.equilDistance();
- const T force = -properDihedral.forceConstant() * properDihedral.multiplicity() * std::sin(deltaPhi);
- const T ePot = properDihedral.forceConstant() * ( 1 + std::cos(deltaPhi) );
+ T deltaPhi = properDihedral.multiplicity() * phi - properDihedral.equilDistance();
+ T force = -properDihedral.forceConstant() * properDihedral.multiplicity() * std::sin(deltaPhi);
+ T ePot = properDihedral.forceConstant() * ( 1 + std::cos(deltaPhi) );
return std::make_tuple(force, ePot);
}
//! \brief Ensure that a geometric quantity lies in (-pi, pi)
-static inline void makeAnglePeriodic(real& angle)
+template<class T>
+inline void makeAnglePeriodic(T& angle)
{
if (angle >= M_PI)
{
}
}
-/*! \brief calculate the cosine of the angle between a and b
+/*! \brief calculate the cosine of the angle between aInput and bInput
*
- * \tparam T float or double
- * \param a a 3D vector
- * \param b another 3D vector
- * \return the cosine of the angle between a and b
+ * \tparam T float or double
+ * \param aInput aInput 3D vector
+ * \param bInput another 3D vector
+ * \return the cosine of the angle between aInput and bInput
*
- * ax*bx + ay*by + az*bz
- * cos-vec (a,b) = ---------------------
- * ||a|| * ||b||
+ * ax*bx + ay*by + az*bz
+ * cos(aInput,bInput) = -----------------------, where aInput = (ax, ay, az)
+ * ||aInput|| * ||bInput||
*/
template<class T>
-inline T basicVectorCosAngle(gmx::BasicVector<T> a, gmx::BasicVector<T> b)
+inline T basicVectorCosAngle(gmx::BasicVector<T> aInput, gmx::BasicVector<T> bInput)
{
- gmx::BasicVector<double> a_double(a[0], a[1], a[2]);
- gmx::BasicVector<double> b_double(b[0], b[1], b[2]);
+ gmx::BasicVector<double> a_double(aInput[0], aInput[1], aInput[2]);
+ gmx::BasicVector<double> b_double(bInput[0], bInput[1], bInput[2]);
double numerator = dot(a_double, b_double);
double denominatorSq = dot(a_double, a_double) * dot(b_double, b_double);
return std::max(cosval, T(-1.0));
}
-//! \brief Computes and returns a dihedral phi angle
-static inline real dihedralPhi(rvec dxIJ, rvec dxKJ, rvec dxKL, rvec m, rvec n)
+/*! \brief compute the angle between vectors a and b
+ *
+ * \tparam T scalar type (float, double, or similar)
+ * \param a a 3D vector
+ * \param b another 3D vector
+ * \return the angle between a and b
+ *
+ * This routine calculates the angle between a & b without any loss of accuracy close to 0/PI.
+ *
+ * Note: This function is not (yet) implemented for the C++ replacement of the
+ * deprecated rvec and dvec.
+ */
+template<class T>
+inline T basicVectorAngle(gmx::BasicVector<T> a, gmx::BasicVector<T> b)
{
- cprod(dxIJ, dxKJ, m);
- cprod(dxKJ, dxKL, n);
- real phi = gmx_angle(m, n);
- real ipr = iprod(dxIJ, n);
- real sign = (ipr < 0.0) ? -1.0 : 1.0;
- phi = sign * phi;
+ gmx::BasicVector<T> w = cross(a, b);
+
+ T wlen = norm(w);
+ T s = dot(a, b);
+
+ return std::atan2(wlen, s);
+}
+
+/*! \brief Computes the dihedral phi angle and two cross products
+ *
+ * \tparam T scalar type (float or double, or similar)
+ * \param[in] dxIJ
+ * \param[in] dxKJ
+ * \param[in] dxKL
+ * \param[out] m output for \p dxIJ x \p dxKJ
+ * \param[out] n output for \p dxKJ x \p dxKL
+ * \return the angle between m and n
+ */
+template<class T>
+inline T dihedralPhi(gmx::BasicVector<T> dxIJ,
+ gmx::BasicVector<T> dxKJ,
+ gmx::BasicVector<T> dxKL,
+ gmx::BasicVector<T>* m,
+ gmx::BasicVector<T>* n)
+{
+ *m = cross(dxIJ, dxKJ);
+ *n = cross(dxKJ, dxKL);
+ T phi = basicVectorAngle(*m, *n);
+ T ipr = dot(dxIJ, *n);
+ T sign = (ipr < 0.0) ? -1.0 : 1.0;
+ phi = sign * phi;
return phi;
}
//! Two-center category common
+//! \brief add shift forces, if requested (static compiler decision)
+template<class T, class ShiftForce>
+inline void
+addShiftForce(const gmx::BasicVector<T>& interactionForce, ShiftForce* shiftForce)
+{
+ *shiftForce += interactionForce;
+}
+
+//! \brief this will be called if shift forces are not computed (and removed by the compiler)
+template<class T>
+inline void addShiftForce([[maybe_unused]] const gmx::BasicVector<T>& fij,
+ [[maybe_unused]] std::nullptr_t*)
+{
+}
+
/*! \brief Spreads and accumulates the bonded forces to the two atoms and adds the virial contribution when needed
*
* \p shiftIndex is used as the periodic shift.
*/
-template <class T>
+template <class T, class ShiftForce>
inline void spreadTwoCenterForces(const T bondForce,
- const gmx::RVec& dx,
- gmx::RVec* force_i,
- gmx::RVec* force_j)
+ const gmx::BasicVector<T>& dx,
+ gmx::BasicVector<T>* force_i,
+ gmx::BasicVector<T>* force_j,
+ ShiftForce* shf_ij,
+ ShiftForce* shf_c)
{
- for (int m = 0; m < dimSize; m++) /* 15 */
- {
- const T fij = bondForce * dx[m];
- (*force_i)[m] += fij;
- (*force_j)[m] -= fij;
- }
+ gmx::BasicVector<T> fij = bondForce * dx;
+ *force_i += fij;
+ *force_j -= fij;
+
+ addShiftForce(fij, shf_ij);
+ addShiftForce(T(-1.0)*fij, shf_c);
+ /* 15 Total */
}
//! Three-center category common
/*! \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
+ * \param cos_theta
+ * \param force
+ * \param r_ij
+ * \param r_kj
+ * \param force_i
+ * \param force_j
+ * \param force_k
*/
-template <class T>
+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,
gmx::BasicVector<T>* force_i,
gmx::BasicVector<T>* force_j,
- gmx::BasicVector<T>* force_k)
+ gmx::BasicVector<T>* force_k,
+ ShiftForce* shf_ij,
+ ShiftForce* shf_kj,
+ ShiftForce* shf_c)
{
T cos_theta2 = cos_theta * cos_theta;
- if (cos_theta2 < 1)
+ if (cos_theta2 < 1) /* 1 */
{
T st = force / std::sqrt(1 - cos_theta2); /* 12 */
- T sth = st * cos_theta; /* 1 */
+ T sth = st * cos_theta; /* 1 */
T nrij2 = dot(r_ij, r_ij); /* 5 */
T nrkj2 = dot(r_kj, r_kj); /* 5 */
- T nrij_1 = 1.0 / std::sqrt(nrij2); /* 10 */
- T nrkj_1 = 1.0 / std::sqrt(nrkj2); /* 10 */
-
- T cik = st * nrij_1 * nrkj_1; /* 2 */
- T cii = sth * nrij_1 * nrij_1; /* 2 */
- T ckk = sth * nrkj_1 * nrkj_1; /* 2 */
-
- gmx::BasicVector<T> f_i{0, 0, 0};
- gmx::BasicVector<T> f_j{0, 0, 0};
- gmx::BasicVector<T> f_k{0, 0, 0};
- for (int m = 0; m < dimSize; m++)
- { /* 39 */
- f_i[m] = -(cik * r_kj[m] - cii * r_ij[m]);
- f_k[m] = -(cik * r_ij[m] - ckk * r_kj[m]);
- f_j[m] = -f_i[m] - f_k[m];
- (*force_i)[m] += f_i[m];
- (*force_j)[m] += f_j[m];
- (*force_k)[m] += f_k[m];
- }
- } /* 161 TOTAL */
+ 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_j = T(-1.0) * (f_i + f_k);
+ *force_i += f_i;
+ *force_j += f_j;
+ *force_k += f_k;
+
+ addShiftForce(f_i, shf_ij);
+ addShiftForce(f_j, shf_c);
+ addShiftForce(f_k, shf_kj);
+
+ } /* 70 TOTAL */
}
//! Four-center category common
-template <class T>
-inline void spreadFourCenterForces(T force, rvec dxIJ, rvec dxJK, rvec dxKL, rvec m, rvec n,
- gmx::RVec* force_i,
- gmx::RVec* force_j,
- gmx::RVec* force_k,
- gmx::RVec* force_l)
+template <class T, class ShiftForce>
+inline void spreadFourCenterForces(T force,
+ const gmx::BasicVector<T>& dxIJ,
+ const gmx::BasicVector<T>& dxJK,
+ const gmx::BasicVector<T>& dxKL,
+ const gmx::BasicVector<T>& m,
+ const gmx::BasicVector<T>& n,
+ gmx::BasicVector<T>* force_i,
+ gmx::BasicVector<T>* force_j,
+ gmx::BasicVector<T>* force_k,
+ gmx::BasicVector<T>* force_l,
+ ShiftForce* shf_ij,
+ ShiftForce* shf_kj,
+ ShiftForce* shf_lj,
+ ShiftForce* shf_c)
{
- rvec f_i, f_j, f_k, f_l;
- rvec uvec, vvec, svec;
- T iprm = iprod(m, m); /* 5 */
- T iprn = iprod(n, n); /* 5 */
- T nrkj2 = iprod(dxJK, dxJK); /* 5 */
- T toler = nrkj2 * GMX_REAL_EPS;
- if ((iprm > toler) && (iprn > toler))
+ T norm2_m = dot(m, m); /* 5 */
+ T norm2_n = dot(n, n); /* 5 */
+ T norm2_jk = dot(dxJK, dxJK); /* 5 */
+ T toler = norm2_jk * GMX_REAL_EPS;
+ if ((norm2_m > toler) && (norm2_n > toler))
{
- T nrkj_1 = gmx::invsqrt(nrkj2); /* 10 */
- T nrkj_2 = nrkj_1 * nrkj_1; /* 1 */
- T nrkj = nrkj2 * nrkj_1; /* 1 */
- T a = -force * nrkj / iprm; /* 11 */
- svmul(a, m, f_i); /* 3 */
- T b = force * nrkj / iprn; /* 11 */
- svmul(b, n, f_l); /* 3 */
- T p = iprod(dxIJ, dxJK); /* 5 */
- p *= nrkj_2; /* 1 */
- T q = iprod(dxKL, dxJK); /* 5 */
- q *= nrkj_2; /* 1 */
- svmul(p, f_i, uvec); /* 3 */
- svmul(q, f_l, vvec); /* 3 */
- rvec_sub(uvec, vvec, svec); /* 3 */
- rvec_sub(f_i, svec, f_j); /* 3 */
- rvec_add(f_l, svec, f_k); /* 3 */
- rvec_inc(force_i->as_vec(), f_i); /* 3 */
- rvec_dec(force_j->as_vec(), f_j); /* 3 */
- rvec_dec(force_k->as_vec(), f_k); /* 3 */
- rvec_inc(force_l->as_vec(), f_l); /* 3 */
+ T rcp_norm2_jk = 1.0f / norm2_jk; /* 1 */
+ T norm_jk = std::sqrt(norm2_jk); /* 10 */
+
+ T a = -force * norm_jk / norm2_m; /* 11 */
+ gmx::BasicVector<T> f_i = a * m; /* 3 */
+
+ T b = force * norm_jk / norm2_n; /* 11 */
+ gmx::BasicVector<T> f_l = b * n; /* 3 */
+
+ T p = rcp_norm2_jk * dot(dxIJ, dxJK); /* 6 */
+ T q = rcp_norm2_jk * dot(dxKL, dxJK); /* 6 */
+ gmx::BasicVector<T> svec = p * f_i - q * f_l; /* 9 */
+
+ gmx::BasicVector<T> f_j = svec - f_i; /* 3 */
+ gmx::BasicVector<T> f_k = T(-1.0)*svec - f_l; /* 6 */
+
+ *force_i += f_i; /* 3 */
+ *force_j += f_j; /* 3 */
+ *force_k += f_k; /* 3 */
+ *force_l += f_l; /* 3 */
+
+ addShiftForce(f_i, shf_ij);
+ addShiftForce(f_j, shf_c);
+ addShiftForce(f_k, shf_kj);
+ addShiftForce(f_l, shf_lj);
}
}
} // namespace nblib
+
+#undef NBLIB_ALWAYS_INLINE
+
#endif // NBLIB_LISTEDFORCES_KERNELS_HPP
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff