From 44faae71bff9b6da158e66c0e14e43c06838cc35 Mon Sep 17 00:00:00 2001 From: Joe Jordan Date: Mon, 16 Aug 2021 14:18:22 +0000 Subject: [PATCH] Updates to nblib listed forces kernels --- api/nblib/listed_forces/dataflow.hpp | 134 +++++++++--- api/nblib/listed_forces/kernels.hpp | 313 +++++++++++++++++---------- 2 files changed, 292 insertions(+), 155 deletions(-) diff --git a/api/nblib/listed_forces/dataflow.hpp b/api/nblib/listed_forces/dataflow.hpp index 1f39f84863..1d20bf704a 100644 --- a/api/nblib/listed_forces/dataflow.hpp +++ b/api/nblib/listed_forces/dataflow.hpp @@ -64,9 +64,24 @@ namespace nblib { -template +//! \brief returns the address of an element in the shiftForces buffer +template +inline ShiftForce* accessShiftForces(int shiftIndex, gmx::ArrayRef 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 shiftForces) +{ + return nullptr; +} + +template 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; @@ -78,7 +93,7 @@ auto computeTwoCenter(const TwoCenterType& parameters, const BasicVector& dx, Ba if (dr2 != 0.0) { force /= dr; - spreadTwoCenterForces(force, dx, fi, fj); + spreadTwoCenterForces(force, dx, fi, fj, sh_f, sh_fc); } return energy; @@ -97,54 +112,59 @@ auto computeTwoCenter(const TwoCenterType& parameters, const BasicVector& dx, Ba * \param[in] pbc Object used for computing distances accounting for PBC's * \return Computed kernel energies */ -template {}>* = nullptr> inline NBLIB_ALWAYS_INLINE auto dispatchInteraction(InteractionIndex index, gmx::ArrayRef bondInstances, gmx::ArrayRef x, Buffer* forces, + gmx::ArrayRef shiftForces, const Pbc& pbc) { KernelEnergy> 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 +template inline NBLIB_ALWAYS_INLINE std::enable_if_t::value, BasicVectorValueType_t> 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 +template inline NBLIB_ALWAYS_INLINE std::enable_if_t::value, BasicVectorValueType_t> addTwoCenterAggregate([[maybe_unused]] const ThreeCenterType& parameters, @@ -152,15 +172,19 @@ 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 +template 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; // calculate 3-center common quantities: angle between x1-x2 and x2-x3 @@ -168,10 +192,10 @@ auto computeThreeCenter(const ThreeCenterType& parameters, const BasicVector& ri 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; } @@ -189,13 +213,14 @@ auto computeThreeCenter(const ThreeCenterType& parameters, const BasicVector& ri * \param[in] PBC * \return Computed kernel energies */ -template {}>* = nullptr> inline NBLIB_ALWAYS_INLINE auto dispatchInteraction(InteractionIndex index, gmx::ArrayRef parameters, gmx::ArrayRef x, Buffer* forces, + gmx::ArrayRef shiftForces, const Pbc& pbc) { KernelEnergy> energy; @@ -212,12 +237,16 @@ auto dispatchInteraction(InteractionIndex index, 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; @@ -276,13 +305,14 @@ addThreeCenterAggregate([[maybe_unused]] const FourCenterType& parameters, * \param[in] pbc Object used for computing distances accounting for PBC's * \return Computed kernel energies */ -template {}>* = nullptr> inline NBLIB_ALWAYS_INLINE auto dispatchInteraction(InteractionIndex index, gmx::ArrayRef parameters, gmx::ArrayRef x, Buffer* forces, + gmx::ArrayRef shiftForces, const Pbc& pbc) { using RealScalar = BasicVectorValueType_t; @@ -300,22 +330,29 @@ auto dispatchInteraction(InteractionIndex index, 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; @@ -338,13 +375,14 @@ auto dispatchInteraction(InteractionIndex index, * \param[in] pbc Object used for computing distances accounting for PBC's * \return Computed kernel energies */ -template {}>* = nullptr> inline NBLIB_ALWAYS_INLINE auto dispatchInteraction(InteractionIndex index, gmx::ArrayRef parameters, gmx::ArrayRef x, Buffer* forces, + [[maybe_unused]] gmx::ArrayRef shiftForces, const Pbc& pbc) { KernelEnergy> energy; @@ -389,23 +427,35 @@ auto dispatchInteraction(InteractionIndex index, * \param[in] pbc Object used for computing distances accounting for PBC's * \return Computed kernel energies */ -template +template auto computeForces(gmx::ArrayRef indices, gmx::ArrayRef parameters, gmx::ArrayRef x, Buffer* forces, + gmx::ArrayRef shiftForces, const Pbc& pbc) { KernelEnergy> 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 +auto computeForces(gmx::ArrayRef indices, + gmx::ArrayRef parameters, + gmx::ArrayRef x, + Buffer* forces, + const Pbc& pbc) +{ + return computeForces(indices, parameters, x, forces, gmx::ArrayRef{}, pbc); +} + /*! \brief implement a loop over bond types and accumulate their force contributions * * \param[in] interactions interaction pairs and bond parameters @@ -414,10 +464,11 @@ auto computeForces(gmx::ArrayRef indices, * \param[in] pbc Object used for computing distances accounting for PBC's * \return Computed kernel energies */ -template +template auto reduceListedForces(const ListedInteractionData& interactions, gmx::ArrayRef x, Buffer* forces, + gmx::ArrayRef shiftForces, const Pbc& pbc) { using ValueType = BasicVectorValueType_t; @@ -425,12 +476,13 @@ auto reduceListedForces(const ListedInteractionData& interactions, 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::type; gmx::ArrayRef> indices(interactionElement.indices); gmx::ArrayRef parameters(interactionElement.parameters); - KernelEnergy energy = computeForces(indices, parameters, x, forces, pbc); + + KernelEnergy energy = computeForces(indices, parameters, x, forces, shiftForces, pbc); energies[CarrierIndex{}] += energy.carrier(); energies[TwoCenterAggregateIndex{}] += energy.twoCenterAggregate(); @@ -444,6 +496,16 @@ auto reduceListedForces(const ListedInteractionData& interactions, return energies; } +//! \brief convenience overload without shift forces +template +auto reduceListedForces(const ListedInteractionData& interactions, + gmx::ArrayRef x, + Buffer* forces, + const Pbc& pbc) +{ + return reduceListedForces(interactions, x, forces, gmx::ArrayRef{}, pbc); +} + } // namespace nblib #undef NBLIB_ALWAYS_INLINE diff --git a/api/nblib/listed_forces/kernels.hpp b/api/nblib/listed_forces/kernels.hpp index ba4b1abfad..b9064ce211 100644 --- a/api/nblib/listed_forces/kernels.hpp +++ b/api/nblib/listed_forces/kernels.hpp @@ -1,7 +1,7 @@ /* * 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. @@ -47,9 +47,11 @@ #include -#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 { @@ -133,12 +135,12 @@ inline std::tuple g96ScalarForce(T k, T x0, T x) 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 @@ -163,8 +165,8 @@ inline std::tuple g96ScalarForce(T kA, T kB, T xA, T xB, T x, T lambda) 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; @@ -173,12 +175,14 @@ inline std::tuple g96ScalarForce(T kA, T kB, T xA, T xB, T x, T lambda) /* 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 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); } @@ -200,11 +204,11 @@ inline std::tuple morseScalarForce(T k, T beta, T x0, T x) 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 @@ -234,7 +238,7 @@ inline std::tuple morseScalarForce(T kA, T kB, T betaA, T betaB, T xA, 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 @@ -271,7 +275,7 @@ inline std::tuple FENEScalarForce(T k, T x0, T x) 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); @@ -284,7 +288,9 @@ inline std::tuple FENEScalarForce(T k, T x0, T x) template 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); } @@ -302,6 +308,7 @@ template inline std::tuple 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; @@ -390,11 +397,8 @@ inline auto bondKernel(T dr, const HalfAttractiveQuarticBondType& bond) //! Three-center interaction type kernels -// linear angles go here -// quartic angles go here - -//! Three-center interaction type dispatch +//! Harmonic Angle template inline auto threeCenterKernel(T dr, const HarmonicAngle& angle) { @@ -406,15 +410,16 @@ inline auto threeCenterKernel(T dr, const HarmonicAngle& angle) template 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 +inline void makeAnglePeriodic(T& angle) { if (angle >= M_PI) { @@ -426,22 +431,22 @@ static inline void makeAnglePeriodic(real& angle) } } -/*! \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 -inline T basicVectorCosAngle(gmx::BasicVector a, gmx::BasicVector b) +inline T basicVectorCosAngle(gmx::BasicVector aInput, gmx::BasicVector bInput) { - gmx::BasicVector a_double(a[0], a[1], a[2]); - gmx::BasicVector b_double(b[0], b[1], b[2]); + gmx::BasicVector a_double(aInput[0], aInput[1], aInput[2]); + gmx::BasicVector 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); @@ -453,15 +458,52 @@ inline T basicVectorCosAngle(gmx::BasicVector a, gmx::BasicVector b) 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 +inline T basicVectorAngle(gmx::BasicVector a, gmx::BasicVector 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 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 +inline T dihedralPhi(gmx::BasicVector dxIJ, + gmx::BasicVector dxKJ, + gmx::BasicVector dxKL, + gmx::BasicVector* m, + gmx::BasicVector* 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; } @@ -501,114 +543,147 @@ inline auto fourCenterKernel(T phi, const RyckaertBellemanDihedral& ryckaertBell //! Two-center category common +//! \brief add shift forces, if requested (static compiler decision) +template +inline void +addShiftForce(const gmx::BasicVector& interactionForce, ShiftForce* shiftForce) +{ + *shiftForce += interactionForce; +} + +//! \brief this will be called if shift forces are not computed (and removed by the compiler) +template +inline void addShiftForce([[maybe_unused]] const gmx::BasicVector& 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 +template inline void spreadTwoCenterForces(const T bondForce, - const gmx::RVec& dx, - gmx::RVec* force_i, - gmx::RVec* force_j) + const gmx::BasicVector& dx, + gmx::BasicVector* force_i, + gmx::BasicVector* 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 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 +template inline void spreadThreeCenterForces(T cos_theta, T force, const gmx::BasicVector& r_ij, const gmx::BasicVector& r_kj, gmx::BasicVector* force_i, gmx::BasicVector* force_j, - gmx::BasicVector* force_k) + gmx::BasicVector* 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 f_i{0, 0, 0}; - gmx::BasicVector f_j{0, 0, 0}; - gmx::BasicVector 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 f_i = cii * r_ij - cik * r_kj; + gmx::BasicVector f_k = ckk * r_kj - cik * r_ij; + gmx::BasicVector 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 -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 +inline void spreadFourCenterForces(T force, + const gmx::BasicVector& dxIJ, + const gmx::BasicVector& dxJK, + const gmx::BasicVector& dxKL, + const gmx::BasicVector& m, + const gmx::BasicVector& n, + gmx::BasicVector* force_i, + gmx::BasicVector* force_j, + gmx::BasicVector* force_k, + gmx::BasicVector* 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 f_i = a * m; /* 3 */ + + T b = force * norm_jk / norm2_n; /* 11 */ + gmx::BasicVector 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 svec = p * f_i - q * f_l; /* 9 */ + + gmx::BasicVector f_j = svec - f_i; /* 3 */ + gmx::BasicVector 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 -- 2.22.0