From b36dd96e4a8ef47460675ddcf9cc904d978b6d59 Mon Sep 17 00:00:00 2001 From: Pascal Merz Date: Wed, 10 Feb 2021 17:27:01 -0700 Subject: [PATCH] Implement velocity terms for virtual sites This change * introduces velocity terms for virtual site reporting replacing the displacement hack used prior, * removes work-arounds used to keep virtual sites frozen during real atom propagation, * allows users of the VirtualSiteHandler object to chose whether to calculate positions only (used for non-output steps), positions and velocities (used for output-steps except in md-vv), or velocities only (used for output-steps in md-vv), which mitigates the fact that calculating the velocities requires a few additional operations, * adds tests of the virtual velocity, * updates the documentation accordingly. Fixes #3866 --- .../functions/interaction-methods.rst | 123 ++- src/gromacs/domdec/domdec.h | 4 +- src/gromacs/domdec/domdec_vsite.cpp | 12 +- src/gromacs/mdlib/update.cpp | 33 +- src/gromacs/mdlib/vsite.cpp | 709 ++++++++++++++---- src/gromacs/mdlib/vsite.h | 19 +- src/gromacs/mdrun/md.cpp | 22 +- src/gromacs/mdrun/mimic.cpp | 2 +- src/gromacs/mdrun/minimize.cpp | 4 +- src/gromacs/mdrun/rerun.cpp | 8 +- src/gromacs/mdrun/shellfc.cpp | 2 +- .../modularsimulator/statepropagatordata.cpp | 4 +- src/programs/mdrun/tests/virtualsites.cpp | 3 +- 13 files changed, 719 insertions(+), 226 deletions(-) diff --git a/docs/reference-manual/functions/interaction-methods.rst b/docs/reference-manual/functions/interaction-methods.rst index 4d88abfb07..524c19f78a 100644 --- a/docs/reference-manual/functions/interaction-methods.rst +++ b/docs/reference-manual/functions/interaction-methods.rst @@ -195,11 +195,16 @@ The force is then redistributed using the same weights: The types of virtual sites supported in |Gromacs| are given in the list below. Constructing atoms in virtual sites can be virtual sites themselves, but only if they are higher in the list, i.e. virtual sites -can be constructed from “particles” that are simpler virtual sites. +can be constructed from “particles” that are simpler virtual sites. The +virtual site velocities are reported, but not used in the integration +of the virtual site positions. -- On top of an atom. This allows giving an atom multiple atom types and +On top of an atom +~~~~~~~~~~~~~~~~~ + +- This allows giving an atom multiple atom types and with that also assigned multiple, different bonded interactions. This - can espically be of use in free-energy calculations. + can especially be of use in free-energy calculations. - The coordinates of the virtual site equal that of the constructing atom: @@ -211,8 +216,15 @@ can be constructed from “particles” that are simpler virtual sites. .. math:: \mathbf{F}_i ~=~ \mathbf{F}_{s} :label: eqnvsite1force -- As a linear combination of two atoms - (:numref:`Fig. %s ` 2): +- The velocity of the virtual site equals that of the constructing atom: + + .. math:: \mathbf{v}_s ~=~ \mathbf{v}_i + :label: eqnvsite1vel + +As a linear combination of two atoms (:numref:`Fig. %s ` 2) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +- The weights are calculated as .. math:: w_i = 1 - a ~,~~ w_j = a :label: eqnvsitelin2atom @@ -220,8 +232,13 @@ can be constructed from “particles” that are simpler virtual sites. - In this case the virtual site is on the line through atoms :math:`i` and :math:`j`. -- On the line through two atoms, with a fixed distance - (:numref:`Fig. %s ` 2fd): +- The velocity of the virtual site is a linear combination of the + velocities of the constructing atoms + +On the line through two atoms, with a fixed distance (:numref:`Fig. %s ` 2fd) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +- The position is calculated as: .. math:: \mathbf{r}_s ~=~ \mathbf{r}_i + a \frac{ \mathbf{r}_{ij} } { | \mathbf{r}_{ij} | } @@ -244,8 +261,18 @@ can be constructed from “particles” that are simpler virtual sites. \end{array} :label: eqnvsite2fdforce -- As a linear combination of three atoms - (:numref:`Fig. %s ` 3): +- The velocity is calculated as: + + .. math:: \mathbf{v}_{s} = \mathbf{v}_{i} + \frac{a}{|\mathbf{r}_{ij}|} + \left(\mathbf{v}_{ij} - \mathbf{r}_{ij} + \frac{\mathbf{v}_{ij}\cdot\mathbf{r}_{ij}} + {|\mathbf{r}_{ij}|^2}\right) + :label: eqnvsite2fdatomvel + +As a linear combination of three atoms (:numref:`Fig. %s ` 3) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +- The weights are calculated as: .. math:: w_i = 1 - a - b ~,~~ w_j = a ~,~~ w_k = b :label: eqnvsitelin3atom @@ -253,11 +280,14 @@ can be constructed from “particles” that are simpler virtual sites. - In this case the virtual site is in the plane of the other three particles. -- In the plane of three atoms, with a fixed distance - (:numref:`Fig. %s ` 3fd): +In the plane of three atoms, with a fixed distance (:numref:`Fig. %s ` 3fd) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - .. math:: \mathbf{r}_s ~=~ \mathbf{r}_i + b \frac{ (1 - a) \mathbf{r}_{ij} + a \mathbf{r}_{jk} } - { | (1 - a) \mathbf{r}_{ij} + a \mathbf{r}_{jk} | } +- The position is calculated as: + + .. math:: \mathbf{r}_s ~=~ \mathbf{r}_i + b \frac{ \mathbf{r}_{ijk} } { | \mathbf{r}_{ijk} | } + ~\mbox{ where }~ + \mathbf{r}_{ijk} ~=~ (1 - a) \mathbf{r}_{ij} + a \mathbf{r}_{jk} :label: eqnvsiteplane3atom - In this case the virtual site is in the plane of the other three @@ -278,8 +308,19 @@ can be constructed from “particles” that are simpler virtual sites. \end{array} :label: eqnvsiteplane3atomforce -- In the plane of three atoms, with a fixed angle and - distance (:numref:`Fig. %s ` 3fad): +- The velocity is calculated as: + + .. math:: \mathbf{v}_{s} ~=~ \mathbf{v}_{i} + + \frac{b}{|\mathbf{r}_{ijk}|} + \left(\dot{\mathbf{r}}_{ijk} - + \mathbf{r}_{ijk}\frac{\dot{\mathbf{r}}_{ijk}\cdot\mathbf{r}_{ijk}} + {|\mathbf{r}_{ijk}|^2}\right) + :label: eqnvsiteplane3atomvel + +In the plane of three atoms, with a fixed angle and distance (:numref:`Fig. %s ` 3fad) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +- The position is calculated as: .. math:: \mathbf{r}_s ~=~ \mathbf{r}_i + d \cos \theta \frac{\mathbf{r}_{ij}}{ | \mathbf{r}_{ij} | } + @@ -334,8 +375,27 @@ can be constructed from “particles” that are simpler virtual sites. \end{array} :label: eqnvsite2fadFforce -- As a non-linear combination of three atoms, out of - plane (:numref:`Fig. %s ` 3out): +- The velocity is calculated as: + + .. math:: \mathbf{v}_{s} &= \mathbf{v}_{i} + d\cos\theta\ \frac{\delta}{\delta t}\frac{\mathbf{r}_{ij}}{|\mathbf{r}_{ij}|} + + d\sin\theta\ \frac{\delta}{\delta t}\frac{\mathbf{r}_{\perp}}{|\mathbf{r}_{\perp}|} \\ + ~\mbox{where}~&\\ + \frac{\delta}{\delta t}\frac{\mathbf{r}_{ij}}{|\mathbf{r}_{ij}|} &= + \frac{1}{|\mathbf{r}_{ij}|}\left(\mathbf{v}_{ij} - \mathbf{r}_{ij} + \frac{\mathbf{v}_{ij}\cdot\mathbf{r}_{ij}}{|\mathbf{r}_{ij}|^2}\right)\\ + \frac{\delta}{\delta t}\frac{\mathbf{r}_{\perp}}{|\mathbf{r}_{\perp}|} &= + \frac{1}{|\mathbf{r}_{\perp}|} + \left(\dot{\mathbf{r}}_{\perp} - \mathbf{r}_{\perp}\frac{\dot{\mathbf{r}}_{\perp}\cdot\mathbf{r}_{\perp}}{|\mathbf{r}_{\perp}|^2}\right)\\ + \dot{\mathbf{r}}_\perp &= \mathbf{v}_{jk} - \mathbf{r}_{ij} + \frac{|\mathbf{r}_{ij}|^2(\mathbf{v}_{ij}\cdot\mathbf{r}_{jk} + \mathbf{r}_{ij}\cdot\mathbf{v}_{jk}) - + (\mathbf{r}_{ij}\cdot\mathbf{r}_{jk})(2\mathbf{r}_{ij}\cdot\mathbf{v}_{ij})} {|\mathbf{r}_{ij}|^4} - + \frac{\mathbf{r}_{ij}\cdot\mathbf{r}_{jk}}{|\mathbf{r}_{ij}|^2}\ \mathbf{v}_{ij} + :label: eqnvsite2fadvel + +As a non-linear combination of three atoms, out of plane (:numref:`Fig. %s ` 3out) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +- The position is calculated as: .. math:: \mathbf{r}_s ~=~ \mathbf{r}_i + a \mathbf{r}_{ij} + b \mathbf{r}_{ik} + c (\mathbf{r}_{ij} \times \mathbf{r}_{ik}) @@ -360,8 +420,15 @@ can be constructed from “particles” that are simpler virtual sites. \end{array} :label: eqnvsitenonlin3atomforce -- From four atoms, with a fixed distance, see - separate :numref:`Fig. %s `. This construction is a bit complex, +- The velocity is calculated as: + + .. math:: \mathbf{v}_{s} ~=~ \mathbf{v}_{i} + \frac{c}{|\mathbf{r}_{m}|}\left(\dot{\mathbf{r}}_{m} - + \mathbf{r}_{m} \frac{\dot{\mathbf{r}}_{m}\cdot\mathbf{r}_{m}}{|\mathbf{r}_{m}|^2}\right) + :label: eqnvsitenonlin3atomvel + +From four atoms, with a fixed distance, see separate :numref:`Fig. %s ` +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +- This construction is a bit complex, in particular since the previous type (4fd) could be unstable which forced us to introduce a more elaborate construction: @@ -373,15 +440,23 @@ can be constructed from “particles” that are simpler virtual sites. The new 4fdn virtual site construction, which is stable even when all constructing atoms are in the same plane. -- +- The position is calculated as + .. math:: \begin{aligned} \mathbf{r}_{ja} &=& a\, \mathbf{r}_{ik} - \mathbf{r}_{ij} = a\, (\mathbf{x}_k - \mathbf{x}_i) - (\mathbf{x}_j - \mathbf{x}_i) \nonumber \\ \mathbf{r}_{jb} &=& b\, \mathbf{r}_{il} - \mathbf{r}_{ij} = b\, (\mathbf{x}_l - \mathbf{x}_i) - (\mathbf{x}_j - \mathbf{x}_i) \nonumber \\ \mathbf{r}_m &=& \mathbf{r}_{ja} \times \mathbf{r}_{jb} \nonumber \\ - \mathbf{x}_s &=& \mathbf{x}_i + c \frac{\mathbf{r}_m}{ | \mathbf{r}_m | } + \mathbf{r}_s &=& \mathbf{r}_i + c \frac{\mathbf{r}_m}{ | \mathbf{r}_m | } \end{aligned} :label: eqnvsite +- The velocity is calculated as: + + .. math:: \mathbf{v}_{s} = \mathbf{v}_{i} + \frac{c}{|\mathbf{r}_{m}|}\left(\dot{\mathbf{r}}_{m} - \mathbf{r}_{m} \frac{\dot{\mathbf{r}}_{m}\cdot\mathbf{r}_{m}}{|\mathbf{r}_{m}|^2}\right)\\ + ~\mbox{where}~&\\ + \dot{\mathbf{r}}_{m} &= \dot{\mathbf{r}}_{ja} \times \mathbf{r}_{jb} + \mathbf{r}_{ja} \times \dot{\mathbf{r}}_{jb} + :label: eqnvsitevel + - In this case the virtual site is at a distance of :math:`|c|` from :math:`i`, while :math:`a` and :math:`b` are parameters. **Note** that the vectors :math:`\mathbf{r}_{ik}` and :math:`\mathbf{r}_{ij}` @@ -400,8 +475,10 @@ can be constructed from “particles” that are simpler virtual sites. value 1), but it should not be used for new simulations. All current |Gromacs| tools will automatically generate type 4fdn instead. -- A linear combination of :math:`N` atoms with relative - weights :math:`a_i`. The weight for atom :math:`i` is: +A linear combination of :math:`N` atoms with relative weights :math:`a_i` +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +- The weight for atom :math:`i` is: .. math:: w_i = a_i \left(\sum_{j=1}^N a_j \right)^{-1} :label: eqnvsiterelweight diff --git a/src/gromacs/domdec/domdec.h b/src/gromacs/domdec/domdec.h index 2b65e1989b..542ced9776 100644 --- a/src/gromacs/domdec/domdec.h +++ b/src/gromacs/domdec/domdec.h @@ -3,7 +3,7 @@ * * Copyright (c) 2005 - 2014, The GROMACS development team. * Copyright (c) 2015,2016,2017,2018,2019 by the GROMACS development team. - * Copyright (c) 2020, by the GROMACS development team, led by + * Copyright (c) 2020,2021, 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. @@ -242,6 +242,8 @@ void dd_move_x_constraints(struct gmx_domdec_t* dd, /*! \brief Communicates the coordinates involved in virtual sites */ void dd_move_x_vsites(const gmx_domdec_t& dd, const matrix box, rvec* x); +/*! \brief Communicates the positions and velocities involved in virtual sites */ +void dd_move_x_and_v_vsites(const gmx_domdec_t& dd, const matrix box, rvec* x, rvec* v); /*! \brief Returns the local atom count array for all constraints * diff --git a/src/gromacs/domdec/domdec_vsite.cpp b/src/gromacs/domdec/domdec_vsite.cpp index cc1880dde5..bde62216ca 100644 --- a/src/gromacs/domdec/domdec_vsite.cpp +++ b/src/gromacs/domdec/domdec_vsite.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 2006,2007,2008,2009,2010 by the GROMACS development team. * Copyright (c) 2012,2013,2014,2015,2016 by the GROMACS development team. - * Copyright (c) 2017,2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2017,2018,2019,2020,2021, 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. @@ -93,6 +93,14 @@ void dd_move_x_vsites(const gmx_domdec_t& dd, const matrix box, rvec* x) } } +void dd_move_x_and_v_vsites(const gmx_domdec_t& dd, const matrix box, rvec* x, rvec* v) +{ + if (dd.vsite_comm) + { + dd_move_x_specat(&dd, dd.vsite_comm, box, x, v, FALSE); + } +} + void dd_clear_local_vsite_indices(gmx_domdec_t* dd) { if (dd->vsite_comm) @@ -143,7 +151,7 @@ int dd_make_local_vsites(gmx_domdec_t* dd, int at_start, gmx::ArrayRefvsite_comm, ga2la_specat, at_start, 1, "vsite", ""); + dd, &ireq, dd->vsite_comm, ga2la_specat, at_start, 2, "vsite", ""); /* Fill in the missing indices */ for (int ftype = 0; ftype < F_NRE; ftype++) diff --git a/src/gromacs/mdlib/update.cpp b/src/gromacs/mdlib/update.cpp index ad61495af8..229a029209 100644 --- a/src/gromacs/mdlib/update.cpp +++ b/src/gromacs/mdlib/update.cpp @@ -255,18 +255,6 @@ void Update::update_temperature_constants(const t_inputrec& inputRecord) return impl_->update_temperature_constants(inputRecord); } -/*! \brief Sets the velocities of virtual sites to zero */ -static void clearVsiteVelocities(int start, int nrend, const unsigned short* particleType, rvec* gmx_restrict v) -{ - for (int a = start; a < nrend; a++) - { - if (particleType[a] == eptVSite) - { - clear_rvec(v[a]); - } - } -} - /*! \brief Sets whether we store the updated velocities */ enum class StoreUpdatedVelocities { @@ -663,19 +651,6 @@ static void do_update_md(int start, } else { - /* Use a simple and thus more efficient integration loop. */ - /* The simple loop does not check for particle type (so it can - * be vectorized), which means we need to clear the velocities - * of virtual sites in advance, when present. Note that vsite - * velocities are computed after update and constraints from - * their displacement. - */ - if (md->haveVsites) - { - /* Note: The overhead of this loop is completely neligible */ - clearVsiteVelocities(start, nrend, md->ptype, v); - } - /* We determine if we have a single T-coupling lambda value for all * atoms. That allows for better SIMD acceleration in the template. * If we do not do temperature coupling (in the run or this step), @@ -811,7 +786,7 @@ static void do_update_vv_vel(int start, for (d = 0; d < DIM; d++) { - if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d]) + if ((ptype[n] != eptShell) && !nFreeze[gf][d]) { v[n][d] = mv1 * (mv1 * v[n][d] + 0.5 * (w_dt * mv2 * f[n][d])); } @@ -862,7 +837,7 @@ static void do_update_vv_pos(int start, for (d = 0; d < DIM; d++) { - if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d]) + if ((ptype[n] != eptShell) && !nFreeze[gf][d]) { xprime[n][d] = mr1 * (mr1 * x[n][d] + mr2 * dt * v[n][d]); } @@ -1045,7 +1020,7 @@ static void doSDUpdateGeneral(const gmx_stochd_t& sd, for (int d = 0; d < DIM; d++) { - if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[freezeGroup][d]) + if ((ptype[n] != eptShell) && !nFreeze[freezeGroup][d]) { if (updateType == SDUpdate::ForcesOnly) { @@ -1185,7 +1160,7 @@ static void do_update_bd(int start, } for (d = 0; (d < DIM); d++) { - if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d]) + if ((ptype[n] != eptShell) && !nFreeze[gf][d]) { if (friction_coefficient != 0) { diff --git a/src/gromacs/mdlib/vsite.cpp b/src/gromacs/mdlib/vsite.cpp index 4a421847e6..28190ef811 100644 --- a/src/gromacs/mdlib/vsite.cpp +++ b/src/gromacs/mdlib/vsite.cpp @@ -4,7 +4,7 @@ * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. * Copyright (c) 2013,2014,2015,2016,2017 The GROMACS development team. - * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020,2021, 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. @@ -242,12 +242,12 @@ public: /*! \brief Create positions of vsite atoms based for the local system * - * \param[in,out] x The coordinates - * \param[in] dt The time step - * \param[in,out] v When != nullptr, velocities for vsites are set as displacement/dt - * \param[in] box The box + * \param[in,out] x The coordinates + * \param[in,out] v The velocities, needed if operation requires it + * \param[in] box The box + * \param[in] operation Whether we calculate positions, velocities, or both */ - void construct(ArrayRef x, real dt, ArrayRef v, const matrix box) const; + void construct(ArrayRef x, ArrayRef v, const matrix box, VSiteOperation operation) const; /*! \brief Spread the force operating on the vsite atoms on the surrounding atoms. * @@ -323,87 +323,167 @@ static inline real inverseNorm(const rvec x) return gmx::invsqrt(iprod(x, x)); } +//! Whether we're calculating the virtual site position +enum class VSiteCalculatePosition +{ + Yes, + No +}; +//! Whether we're calculating the virtual site velocity +enum class VSiteCalculateVelocity +{ + Yes, + No +}; + #ifndef DOXYGEN /* Vsite construction routines */ -static void constr_vsite1(const rvec xi, rvec x) -{ - copy_rvec(xi, x); +// GCC 8 falsely flags unused variables if constexpr prunes a code path, fixed in GCC 9 +// https://gcc.gnu.org/bugzilla/show_bug.cgi?id=85827 +// clang-format off +GCC_DIAGNOSTIC_IGNORE(-Wunused-but-set-parameter) +// clang-format on - /* TOTAL: 0 flops */ +template +static void constr_vsite1(const rvec xi, rvec x, const rvec vi, rvec v) +{ + if (calculatePosition == VSiteCalculatePosition::Yes) + { + copy_rvec(xi, x); + /* TOTAL: 0 flops */ + } + if (calculateVelocity == VSiteCalculateVelocity::Yes) + { + copy_rvec(vi, v); + } } -static void constr_vsite2(const rvec xi, const rvec xj, rvec x, real a, const t_pbc* pbc) +template +static void +constr_vsite2(const rvec xi, const rvec xj, rvec x, real a, const t_pbc* pbc, const rvec vi, const rvec vj, rvec v) { - real b = 1 - a; + const real b = 1 - a; /* 1 flop */ - if (pbc) + if (calculatePosition == VSiteCalculatePosition::Yes) { - rvec dx; - pbc_dx_aiuc(pbc, xj, xi, dx); - x[XX] = xi[XX] + a * dx[XX]; - x[YY] = xi[YY] + a * dx[YY]; - x[ZZ] = xi[ZZ] + a * dx[ZZ]; + if (pbc) + { + rvec dx; + pbc_dx_aiuc(pbc, xj, xi, dx); + x[XX] = xi[XX] + a * dx[XX]; + x[YY] = xi[YY] + a * dx[YY]; + x[ZZ] = xi[ZZ] + a * dx[ZZ]; + } + else + { + x[XX] = b * xi[XX] + a * xj[XX]; + x[YY] = b * xi[YY] + a * xj[YY]; + x[ZZ] = b * xi[ZZ] + a * xj[ZZ]; + /* 9 Flops */ + } + /* TOTAL: 10 flops */ } - else + if (calculateVelocity == VSiteCalculateVelocity::Yes) { - x[XX] = b * xi[XX] + a * xj[XX]; - x[YY] = b * xi[YY] + a * xj[YY]; - x[ZZ] = b * xi[ZZ] + a * xj[ZZ]; - /* 9 Flops */ + v[XX] = b * vi[XX] + a * vj[XX]; + v[YY] = b * vi[YY] + a * vj[YY]; + v[ZZ] = b * vi[ZZ] + a * vj[ZZ]; } - - /* TOTAL: 10 flops */ } -static void constr_vsite2FD(const rvec xi, const rvec xj, rvec x, real a, const t_pbc* pbc) +template +static void +constr_vsite2FD(const rvec xi, const rvec xj, rvec x, real a, const t_pbc* pbc, const rvec vi, const rvec vj, rvec v) { - rvec xij; + rvec xij = { 0 }; pbc_rvec_sub(pbc, xj, xi, xij); /* 3 flops */ - const real b = a * inverseNorm(xij); + const real invNormXij = inverseNorm(xij); + const real b = a * invNormXij; /* 6 + 10 flops */ - x[XX] = xi[XX] + b * xij[XX]; - x[YY] = xi[YY] + b * xij[YY]; - x[ZZ] = xi[ZZ] + b * xij[ZZ]; - /* 6 Flops */ + if (calculatePosition == VSiteCalculatePosition::Yes) + { + x[XX] = xi[XX] + b * xij[XX]; + x[YY] = xi[YY] + b * xij[YY]; + x[ZZ] = xi[ZZ] + b * xij[ZZ]; + /* 6 Flops */ + /* TOTAL: 25 flops */ + } + if (calculateVelocity == VSiteCalculateVelocity::Yes) + { + rvec vij = { 0 }; + rvec_sub(vj, vi, vij); + const real vijDotXij = iprod(vij, xij); - /* TOTAL: 25 flops */ + v[XX] = vi[XX] + b * (vij[XX] - xij[XX] * vijDotXij * invNormXij * invNormXij); + v[YY] = vi[YY] + b * (vij[YY] - xij[YY] * vijDotXij * invNormXij * invNormXij); + v[ZZ] = vi[ZZ] + b * (vij[ZZ] - xij[ZZ] * vijDotXij * invNormXij * invNormXij); + } } -static void constr_vsite3(const rvec xi, const rvec xj, const rvec xk, rvec x, real a, real b, const t_pbc* pbc) +template +static void constr_vsite3(const rvec xi, + const rvec xj, + const rvec xk, + rvec x, + real a, + real b, + const t_pbc* pbc, + const rvec vi, + const rvec vj, + const rvec vk, + rvec v) { - real c = 1 - a - b; + const real c = 1 - a - b; /* 2 flops */ - if (pbc) + if (calculatePosition == VSiteCalculatePosition::Yes) { - rvec dxj, dxk; + if (pbc) + { + rvec dxj, dxk; - pbc_dx_aiuc(pbc, xj, xi, dxj); - pbc_dx_aiuc(pbc, xk, xi, dxk); - x[XX] = xi[XX] + a * dxj[XX] + b * dxk[XX]; - x[YY] = xi[YY] + a * dxj[YY] + b * dxk[YY]; - x[ZZ] = xi[ZZ] + a * dxj[ZZ] + b * dxk[ZZ]; + pbc_dx_aiuc(pbc, xj, xi, dxj); + pbc_dx_aiuc(pbc, xk, xi, dxk); + x[XX] = xi[XX] + a * dxj[XX] + b * dxk[XX]; + x[YY] = xi[YY] + a * dxj[YY] + b * dxk[YY]; + x[ZZ] = xi[ZZ] + a * dxj[ZZ] + b * dxk[ZZ]; + } + else + { + x[XX] = c * xi[XX] + a * xj[XX] + b * xk[XX]; + x[YY] = c * xi[YY] + a * xj[YY] + b * xk[YY]; + x[ZZ] = c * xi[ZZ] + a * xj[ZZ] + b * xk[ZZ]; + /* 15 Flops */ + } + /* TOTAL: 17 flops */ } - else + if (calculateVelocity == VSiteCalculateVelocity::Yes) { - x[XX] = c * xi[XX] + a * xj[XX] + b * xk[XX]; - x[YY] = c * xi[YY] + a * xj[YY] + b * xk[YY]; - x[ZZ] = c * xi[ZZ] + a * xj[ZZ] + b * xk[ZZ]; - /* 15 Flops */ + v[XX] = c * vi[XX] + a * vj[XX] + b * vk[XX]; + v[YY] = c * vi[YY] + a * vj[YY] + b * vk[YY]; + v[ZZ] = c * vi[ZZ] + a * vj[ZZ] + b * vk[ZZ]; } - - /* TOTAL: 17 flops */ } -static void constr_vsite3FD(const rvec xi, const rvec xj, const rvec xk, rvec x, real a, real b, const t_pbc* pbc) +template +static void constr_vsite3FD(const rvec xi, + const rvec xj, + const rvec xk, + rvec x, + real a, + real b, + const t_pbc* pbc, + const rvec vi, + const rvec vj, + const rvec vk, + rvec v) { rvec xij, xjk, temp; - real c; pbc_rvec_sub(pbc, xj, xi, xij); pbc_rvec_sub(pbc, xk, xj, xjk); @@ -415,45 +495,122 @@ static void constr_vsite3FD(const rvec xi, const rvec xj, const rvec xk, rvec x, temp[ZZ] = xij[ZZ] + a * xjk[ZZ]; /* 6 flops */ - c = b * inverseNorm(temp); + const real invNormTemp = inverseNorm(temp); + const real c = b * invNormTemp; /* 6 + 10 flops */ - x[XX] = xi[XX] + c * temp[XX]; - x[YY] = xi[YY] + c * temp[YY]; - x[ZZ] = xi[ZZ] + c * temp[ZZ]; - /* 6 Flops */ + if (calculatePosition == VSiteCalculatePosition::Yes) + { + x[XX] = xi[XX] + c * temp[XX]; + x[YY] = xi[YY] + c * temp[YY]; + x[ZZ] = xi[ZZ] + c * temp[ZZ]; + /* 6 Flops */ + /* TOTAL: 34 flops */ + } + if (calculateVelocity == VSiteCalculateVelocity::Yes) + { + rvec vij = { 0 }; + rvec vjk = { 0 }; + rvec_sub(vj, vi, vij); + rvec_sub(vk, vj, vjk); + const rvec tempV = { vij[XX] + a * vjk[XX], vij[YY] + a * vjk[YY], vij[ZZ] + a * vjk[ZZ] }; + const real tempDotTempV = iprod(temp, tempV); - /* TOTAL: 34 flops */ + v[XX] = vi[XX] + c * (tempV[XX] - temp[XX] * tempDotTempV * invNormTemp * invNormTemp); + v[YY] = vi[YY] + c * (tempV[YY] - temp[YY] * tempDotTempV * invNormTemp * invNormTemp); + v[ZZ] = vi[ZZ] + c * (tempV[ZZ] - temp[ZZ] * tempDotTempV * invNormTemp * invNormTemp); + } } -static void constr_vsite3FAD(const rvec xi, const rvec xj, const rvec xk, rvec x, real a, real b, const t_pbc* pbc) -{ +template +static void constr_vsite3FAD(const rvec xi, + const rvec xj, + const rvec xk, + rvec x, + real a, + real b, + const t_pbc* pbc, + const rvec vi, + const rvec vj, + const rvec vk, + rvec v) +{ // Note: a = d * cos(theta) + // b = d * sin(theta) rvec xij, xjk, xp; - real a1, b1, c1, invdij; pbc_rvec_sub(pbc, xj, xi, xij); pbc_rvec_sub(pbc, xk, xj, xjk); /* 6 flops */ - invdij = inverseNorm(xij); - c1 = invdij * invdij * iprod(xij, xjk); - xp[XX] = xjk[XX] - c1 * xij[XX]; - xp[YY] = xjk[YY] - c1 * xij[YY]; - xp[ZZ] = xjk[ZZ] - c1 * xij[ZZ]; - a1 = a * invdij; - b1 = b * inverseNorm(xp); + const real invdij = inverseNorm(xij); + const real xijDotXjk = iprod(xij, xjk); + const real c1 = invdij * invdij * xijDotXjk; + xp[XX] = xjk[XX] - c1 * xij[XX]; + xp[YY] = xjk[YY] - c1 * xij[YY]; + xp[ZZ] = xjk[ZZ] - c1 * xij[ZZ]; + const real a1 = a * invdij; + const real invNormXp = inverseNorm(xp); + const real b1 = b * invNormXp; /* 45 */ - x[XX] = xi[XX] + a1 * xij[XX] + b1 * xp[XX]; - x[YY] = xi[YY] + a1 * xij[YY] + b1 * xp[YY]; - x[ZZ] = xi[ZZ] + a1 * xij[ZZ] + b1 * xp[ZZ]; - /* 12 Flops */ - - /* TOTAL: 63 flops */ + if (calculatePosition == VSiteCalculatePosition::Yes) + { + x[XX] = xi[XX] + a1 * xij[XX] + b1 * xp[XX]; + x[YY] = xi[YY] + a1 * xij[YY] + b1 * xp[YY]; + x[ZZ] = xi[ZZ] + a1 * xij[ZZ] + b1 * xp[ZZ]; + /* 12 Flops */ + /* TOTAL: 63 flops */ + } + + if (calculateVelocity == VSiteCalculateVelocity::Yes) + { + rvec vij = { 0 }; + rvec vjk = { 0 }; + rvec_sub(vj, vi, vij); + rvec_sub(vk, vj, vjk); + + const real vijDotXjkPlusXijDotVjk = iprod(vij, xjk) + iprod(xij, vjk); + const real xijDotVij = iprod(xij, vij); + const real invNormXij2 = invdij * invdij; + + rvec vp = { 0 }; + vp[XX] = vjk[XX] + - xij[XX] * invNormXij2 + * (vijDotXjkPlusXijDotVjk - invNormXij2 * xijDotXjk * xijDotVij * 2) + - vij[XX] * xijDotXjk * invNormXij2; + vp[YY] = vjk[YY] + - xij[YY] * invNormXij2 + * (vijDotXjkPlusXijDotVjk - invNormXij2 * xijDotXjk * xijDotVij * 2) + - vij[YY] * xijDotXjk * invNormXij2; + vp[ZZ] = vjk[ZZ] + - xij[ZZ] * invNormXij2 + * (vijDotXjkPlusXijDotVjk - invNormXij2 * xijDotXjk * xijDotVij * 2) + - vij[ZZ] * xijDotXjk * invNormXij2; + + const real xpDotVp = iprod(xp, vp); + + v[XX] = vi[XX] + a1 * (vij[XX] - xij[XX] * xijDotVij * invdij * invdij) + + b1 * (vp[XX] - xp[XX] * xpDotVp * invNormXp * invNormXp); + v[YY] = vi[YY] + a1 * (vij[YY] - xij[YY] * xijDotVij * invdij * invdij) + + b1 * (vp[YY] - xp[YY] * xpDotVp * invNormXp * invNormXp); + v[ZZ] = vi[ZZ] + a1 * (vij[ZZ] - xij[ZZ] * xijDotVij * invdij * invdij) + + b1 * (vp[ZZ] - xp[ZZ] * xpDotVp * invNormXp * invNormXp); + } } -static void -constr_vsite3OUT(const rvec xi, const rvec xj, const rvec xk, rvec x, real a, real b, real c, const t_pbc* pbc) +template +static void constr_vsite3OUT(const rvec xi, + const rvec xj, + const rvec xk, + rvec x, + real a, + real b, + real c, + const t_pbc* pbc, + const rvec vi, + const rvec vj, + const rvec vk, + rvec v) { rvec xij, xik, temp; @@ -462,14 +619,34 @@ constr_vsite3OUT(const rvec xi, const rvec xj, const rvec xk, rvec x, real a, re cprod(xij, xik, temp); /* 15 Flops */ - x[XX] = xi[XX] + a * xij[XX] + b * xik[XX] + c * temp[XX]; - x[YY] = xi[YY] + a * xij[YY] + b * xik[YY] + c * temp[YY]; - x[ZZ] = xi[ZZ] + a * xij[ZZ] + b * xik[ZZ] + c * temp[ZZ]; - /* 18 Flops */ + if (calculatePosition == VSiteCalculatePosition::Yes) + { + x[XX] = xi[XX] + a * xij[XX] + b * xik[XX] + c * temp[XX]; + x[YY] = xi[YY] + a * xij[YY] + b * xik[YY] + c * temp[YY]; + x[ZZ] = xi[ZZ] + a * xij[ZZ] + b * xik[ZZ] + c * temp[ZZ]; + /* 18 Flops */ + /* TOTAL: 33 flops */ + } + + if (calculateVelocity == VSiteCalculateVelocity::Yes) + { + rvec vij = { 0 }; + rvec vik = { 0 }; + rvec_sub(vj, vi, vij); + rvec_sub(vk, vi, vik); + + rvec temp1 = { 0 }; + rvec temp2 = { 0 }; + cprod(vij, xik, temp1); + cprod(xij, vik, temp2); - /* TOTAL: 33 flops */ + v[XX] = vi[XX] + a * vij[XX] + b * vik[XX] + c * (temp1[XX] + temp2[XX]); + v[YY] = vi[YY] + a * vij[YY] + b * vik[YY] + c * (temp1[YY] + temp2[YY]); + v[ZZ] = vi[ZZ] + a * vij[ZZ] + b * vik[ZZ] + c * (temp1[ZZ] + temp2[ZZ]); + } } +template static void constr_vsite4FD(const rvec xi, const rvec xj, const rvec xk, @@ -478,7 +655,12 @@ static void constr_vsite4FD(const rvec xi, real a, real b, real c, - const t_pbc* pbc) + const t_pbc* pbc, + const rvec vi, + const rvec vj, + const rvec vk, + const rvec vl, + rvec v) { rvec xij, xjk, xjl, temp; real d; @@ -494,17 +676,41 @@ static void constr_vsite4FD(const rvec xi, temp[ZZ] = xij[ZZ] + a * xjk[ZZ] + b * xjl[ZZ]; /* 12 flops */ - d = c * inverseNorm(temp); + const real invRm = inverseNorm(temp); + d = c * invRm; /* 6 + 10 flops */ - x[XX] = xi[XX] + d * temp[XX]; - x[YY] = xi[YY] + d * temp[YY]; - x[ZZ] = xi[ZZ] + d * temp[ZZ]; - /* 6 Flops */ + if (calculatePosition == VSiteCalculatePosition::Yes) + { + x[XX] = xi[XX] + d * temp[XX]; + x[YY] = xi[YY] + d * temp[YY]; + x[ZZ] = xi[ZZ] + d * temp[ZZ]; + /* 6 Flops */ + /* TOTAL: 43 flops */ + } + if (calculateVelocity == VSiteCalculateVelocity::Yes) + { + rvec vij = { 0 }; + rvec vjk = { 0 }; + rvec vjl = { 0 }; - /* TOTAL: 43 flops */ + rvec_sub(vj, vi, vij); + rvec_sub(vk, vj, vjk); + rvec_sub(vl, vj, vjl); + + rvec vm = { 0 }; + vm[XX] = vij[XX] + a * vjk[XX] + b * vjl[XX]; + vm[YY] = vij[YY] + a * vjk[YY] + b * vjl[YY]; + vm[ZZ] = vij[ZZ] + a * vjk[ZZ] + b * vjl[ZZ]; + + const real vmDotRm = iprod(vm, temp); + v[XX] = vi[XX] + d * (vm[XX] - temp[XX] * vmDotRm * invRm * invRm); + v[YY] = vi[YY] + d * (vm[YY] - temp[YY] * vmDotRm * invRm * invRm); + v[ZZ] = vi[ZZ] + d * (vm[ZZ] - temp[ZZ] * vmDotRm * invRm * invRm); + } } +template static void constr_vsite4FDN(const rvec xi, const rvec xj, const rvec xk, @@ -513,7 +719,12 @@ static void constr_vsite4FDN(const rvec xi, real a, real b, real c, - const t_pbc* pbc) + const t_pbc* pbc, + const rvec vi, + const rvec vj, + const rvec vk, + const rvec vl, + rvec v) { rvec xij, xik, xil, ra, rb, rja, rjb, rm; real d; @@ -540,53 +751,121 @@ static void constr_vsite4FDN(const rvec xi, cprod(rja, rjb, rm); /* 9 flops */ - d = c * inverseNorm(rm); + const real invNormRm = inverseNorm(rm); + d = c * invNormRm; /* 5+5+1 flops */ - x[XX] = xi[XX] + d * rm[XX]; - x[YY] = xi[YY] + d * rm[YY]; - x[ZZ] = xi[ZZ] + d * rm[ZZ]; - /* 6 Flops */ + if (calculatePosition == VSiteCalculatePosition::Yes) + { + x[XX] = xi[XX] + d * rm[XX]; + x[YY] = xi[YY] + d * rm[YY]; + x[ZZ] = xi[ZZ] + d * rm[ZZ]; + /* 6 Flops */ + /* TOTAL: 47 flops */ + } - /* TOTAL: 47 flops */ + if (calculateVelocity == VSiteCalculateVelocity::Yes) + { + rvec vij = { 0 }; + rvec vik = { 0 }; + rvec vil = { 0 }; + rvec_sub(vj, vi, vij); + rvec_sub(vk, vi, vik); + rvec_sub(vl, vi, vil); + + rvec vja = { 0 }; + rvec vjb = { 0 }; + + vja[XX] = a * vik[XX] - vij[XX]; + vja[YY] = a * vik[YY] - vij[YY]; + vja[ZZ] = a * vik[ZZ] - vij[ZZ]; + vjb[XX] = b * vil[XX] - vij[XX]; + vjb[YY] = b * vil[YY] - vij[YY]; + vjb[ZZ] = b * vil[ZZ] - vij[ZZ]; + + rvec temp1 = { 0 }; + rvec temp2 = { 0 }; + cprod(vja, rjb, temp1); + cprod(rja, vjb, temp2); + + rvec vm = { 0 }; + vm[XX] = temp1[XX] + temp2[XX]; + vm[YY] = temp1[YY] + temp2[YY]; + vm[ZZ] = temp1[ZZ] + temp2[ZZ]; + + const real rmDotVm = iprod(rm, vm); + v[XX] = vi[XX] + d * (vm[XX] - rm[XX] * rmDotVm * invNormRm * invNormRm); + v[YY] = vi[YY] + d * (vm[YY] - rm[YY] * rmDotVm * invNormRm * invNormRm); + v[ZZ] = vi[ZZ] + d * (vm[ZZ] - rm[ZZ] * rmDotVm * invNormRm * invNormRm); + } } -static int constr_vsiten(const t_iatom* ia, ArrayRef ip, ArrayRef x, const t_pbc* pbc) +template +static int constr_vsiten(const t_iatom* ia, + ArrayRef ip, + ArrayRef x, + const t_pbc* pbc, + ArrayRef v) { rvec x1, dx; dvec dsum; - int n3, av, ai; real a; + dvec dvsum = { 0 }; + rvec v1 = { 0 }; - n3 = 3 * ip[ia[0]].vsiten.n; - av = ia[1]; - ai = ia[2]; + const int n3 = 3 * ip[ia[0]].vsiten.n; + const int av = ia[1]; + int ai = ia[2]; copy_rvec(x[ai], x1); + copy_rvec(v[ai], v1); clear_dvec(dsum); for (int i = 3; i < n3; i += 3) { ai = ia[i + 2]; a = ip[ia[i]].vsiten.a; - if (pbc) + if (calculatePosition == VSiteCalculatePosition::Yes) { - pbc_dx_aiuc(pbc, x[ai], x1, dx); + if (pbc) + { + pbc_dx_aiuc(pbc, x[ai], x1, dx); + } + else + { + rvec_sub(x[ai], x1, dx); + } + dsum[XX] += a * dx[XX]; + dsum[YY] += a * dx[YY]; + dsum[ZZ] += a * dx[ZZ]; + /* 9 Flops */ } - else + if (calculateVelocity == VSiteCalculateVelocity::Yes) { - rvec_sub(x[ai], x1, dx); + rvec_sub(v[ai], v1, dx); + dvsum[XX] += a * dx[XX]; + dvsum[YY] += a * dx[YY]; + dvsum[ZZ] += a * dx[ZZ]; + /* 9 Flops */ } - dsum[XX] += a * dx[XX]; - dsum[YY] += a * dx[YY]; - dsum[ZZ] += a * dx[ZZ]; - /* 9 Flops */ } - x[av][XX] = x1[XX] + dsum[XX]; - x[av][YY] = x1[YY] + dsum[YY]; - x[av][ZZ] = x1[ZZ] + dsum[ZZ]; + if (calculatePosition == VSiteCalculatePosition::Yes) + { + x[av][XX] = x1[XX] + dsum[XX]; + x[av][YY] = x1[YY] + dsum[YY]; + x[av][ZZ] = x1[ZZ] + dsum[ZZ]; + } + + if (calculateVelocity == VSiteCalculateVelocity::Yes) + { + v[av][XX] = v1[XX] + dvsum[XX]; + v[av][YY] = v1[YY] + dvsum[YY]; + v[av][ZZ] = v1[ZZ] + dvsum[ZZ]; + } return n3; } +// End GCC 8 bug +GCC_DIAGNOSTIC_RESET #endif // DOXYGEN @@ -615,29 +894,49 @@ static PbcMode getPbcMode(const t_pbc* pbcPtr) /*! \brief Executes the vsite construction task for a single thread * + * \tparam operation Whether we are calculating positions, velocities, or both * \param[in,out] x Coordinates to construct vsites for - * \param[in] dt Time step, needed when v is not empty - * \param[in,out] v When not empty, velocities are generated for virtual sites + * \param[in,out] v Velocities are generated for virtual sites if `operation` requires it * \param[in] ip Interaction parameters for all interaction, only vsite parameters are used * \param[in] ilist The interaction lists, only vsites are usesd * \param[in] pbc_null PBC struct, used for PBC distance calculations when !=nullptr */ +template static void construct_vsites_thread(ArrayRef x, - const real dt, ArrayRef v, ArrayRef ip, ArrayRef ilist, const t_pbc* pbc_null) { - real inv_dt; - if (!v.empty()) - { - inv_dt = 1.0 / dt; - } - else - { - inv_dt = 1.0; - } + if (calculateVelocity == VSiteCalculateVelocity::Yes) + { + GMX_RELEASE_ASSERT(!v.empty(), + "Can't calculate velocities without access to velocity vector."); + } + + // Work around clang bug (unfixed as of Feb 2021) + // https://bugs.llvm.org/show_bug.cgi?id=35450 + // clang-format off + CLANG_DIAGNOSTIC_IGNORE(-Wunused-lambda-capture) + // clang-format on + // GCC 8 falsely flags unused variables if constexpr prunes a code path, fixed in GCC 9 + // https://gcc.gnu.org/bugzilla/show_bug.cgi?id=85827 + // clang-format off + GCC_DIAGNOSTIC_IGNORE(-Wunused-but-set-parameter) + // clang-format on + // getVOrNull returns a velocity rvec if we need it, nullptr otherwise. + auto getVOrNull = [v](int idx) -> real* { + if (calculateVelocity == VSiteCalculateVelocity::Yes) + { + return v[idx].as_vec(); + } + else + { + return nullptr; + } + }; + GCC_DIAGNOSTIC_RESET + CLANG_DIAGNOSTIC_RESET const PbcMode pbcMode = getPbcMode(pbc_null); /* We need another pbc pointer, as with charge groups we switch per vsite */ @@ -674,39 +973,97 @@ static void construct_vsites_thread(ArrayRef x, real b1, c1; switch (ftype) { - case F_VSITE1: constr_vsite1(x[ai], x[avsite]); break; + case F_VSITE1: + constr_vsite1( + x[ai], x[avsite], getVOrNull(ai), getVOrNull(avsite)); + break; case F_VSITE2: aj = ia[3]; - constr_vsite2(x[ai], x[aj], x[avsite], a1, pbc_null2); + constr_vsite2(x[ai], + x[aj], + x[avsite], + a1, + pbc_null2, + getVOrNull(ai), + getVOrNull(aj), + getVOrNull(avsite)); break; case F_VSITE2FD: aj = ia[3]; - constr_vsite2FD(x[ai], x[aj], x[avsite], a1, pbc_null2); + constr_vsite2FD(x[ai], + x[aj], + x[avsite], + a1, + pbc_null2, + getVOrNull(ai), + getVOrNull(aj), + getVOrNull(avsite)); break; case F_VSITE3: aj = ia[3]; ak = ia[4]; b1 = ip[tp].vsite.b; - constr_vsite3(x[ai], x[aj], x[ak], x[avsite], a1, b1, pbc_null2); + constr_vsite3(x[ai], + x[aj], + x[ak], + x[avsite], + a1, + b1, + pbc_null2, + getVOrNull(ai), + getVOrNull(aj), + getVOrNull(ak), + getVOrNull(avsite)); break; case F_VSITE3FD: aj = ia[3]; ak = ia[4]; b1 = ip[tp].vsite.b; - constr_vsite3FD(x[ai], x[aj], x[ak], x[avsite], a1, b1, pbc_null2); + constr_vsite3FD(x[ai], + x[aj], + x[ak], + x[avsite], + a1, + b1, + pbc_null2, + getVOrNull(ai), + getVOrNull(aj), + getVOrNull(ak), + getVOrNull(avsite)); break; case F_VSITE3FAD: aj = ia[3]; ak = ia[4]; b1 = ip[tp].vsite.b; - constr_vsite3FAD(x[ai], x[aj], x[ak], x[avsite], a1, b1, pbc_null2); + constr_vsite3FAD(x[ai], + x[aj], + x[ak], + x[avsite], + a1, + b1, + pbc_null2, + getVOrNull(ai), + getVOrNull(aj), + getVOrNull(ak), + getVOrNull(avsite)); break; case F_VSITE3OUT: aj = ia[3]; ak = ia[4]; b1 = ip[tp].vsite.b; c1 = ip[tp].vsite.c; - constr_vsite3OUT(x[ai], x[aj], x[ak], x[avsite], a1, b1, c1, pbc_null2); + constr_vsite3OUT(x[ai], + x[aj], + x[ak], + x[avsite], + a1, + b1, + c1, + pbc_null2, + getVOrNull(ai), + getVOrNull(aj), + getVOrNull(ak), + getVOrNull(avsite)); break; case F_VSITE4FD: aj = ia[3]; @@ -714,7 +1071,20 @@ static void construct_vsites_thread(ArrayRef x, al = ia[5]; b1 = ip[tp].vsite.b; c1 = ip[tp].vsite.c; - constr_vsite4FD(x[ai], x[aj], x[ak], x[al], x[avsite], a1, b1, c1, pbc_null2); + constr_vsite4FD(x[ai], + x[aj], + x[ak], + x[al], + x[avsite], + a1, + b1, + c1, + pbc_null2, + getVOrNull(ai), + getVOrNull(aj), + getVOrNull(ak), + getVOrNull(al), + getVOrNull(avsite)); break; case F_VSITE4FDN: aj = ia[3]; @@ -722,9 +1092,24 @@ static void construct_vsites_thread(ArrayRef x, al = ia[5]; b1 = ip[tp].vsite.b; c1 = ip[tp].vsite.c; - constr_vsite4FDN(x[ai], x[aj], x[ak], x[al], x[avsite], a1, b1, c1, pbc_null2); + constr_vsite4FDN(x[ai], + x[aj], + x[ak], + x[al], + x[avsite], + a1, + b1, + c1, + pbc_null2, + getVOrNull(ai), + getVOrNull(aj), + getVOrNull(ak), + getVOrNull(al), + getVOrNull(avsite)); + break; + case F_VSITEN: + inc = constr_vsiten(ia, ip, x, pbc_null2, v); break; - case F_VSITEN: inc = constr_vsiten(ia, ip, x, pbc_null2); break; default: gmx_fatal(FARGS, "No such vsite type %d in %s, line %d", ftype, __FILE__, __LINE__); } @@ -739,13 +1124,6 @@ static void construct_vsites_thread(ArrayRef x, rvec_add(xv, dx, x[avsite]); } } - if (!v.empty()) - { - /* Calculate velocity of vsite... */ - rvec vv; - rvec_sub(x[avsite], xv, vv); - svmul(inv_dt, vv, v[avsite]); - } /* Increment loop variables */ i += inc; @@ -759,16 +1137,15 @@ static void construct_vsites_thread(ArrayRef x, * * \param[in] threadingInfo Used to divide work over threads when != nullptr * \param[in,out] x Coordinates to construct vsites for - * \param[in] dt Time step, needed when v is not empty * \param[in,out] v When not empty, velocities are generated for virtual sites * \param[in] ip Interaction parameters for all interaction, only vsite parameters are used * \param[in] ilist The interaction lists, only vsites are usesd * \param[in] domainInfo Information about PBC and DD * \param[in] box Used for PBC when PBC is set in domainInfo */ +template static void construct_vsites(const ThreadingInfo* threadingInfo, ArrayRef x, - real dt, ArrayRef v, ArrayRef ip, ArrayRef ilist, @@ -801,12 +1178,20 @@ static void construct_vsites(const ThreadingInfo* threadingInfo, if (useDomdec) { - dd_move_x_vsites(*domainInfo.domdec_, box, as_rvec_array(x.data())); + if (calculateVelocity == VSiteCalculateVelocity::Yes) + { + dd_move_x_and_v_vsites( + *domainInfo.domdec_, box, as_rvec_array(x.data()), as_rvec_array(v.data())); + } + else + { + dd_move_x_vsites(*domainInfo.domdec_, box, as_rvec_array(x.data())); + } } if (threadingInfo == nullptr || threadingInfo->numThreads() == 1) { - construct_vsites_thread(x, dt, v, ip, ilist, pbc_null); + construct_vsites_thread(x, v, ip, ilist, pbc_null); } else { @@ -819,31 +1204,52 @@ static void construct_vsites(const ThreadingInfo* threadingInfo, GMX_ASSERT(tData.rangeStart >= 0, "The thread data should be initialized before calling construct_vsites"); - construct_vsites_thread(x, dt, v, ip, tData.ilist, pbc_null); + construct_vsites_thread( + x, v, ip, tData.ilist, pbc_null); if (tData.useInterdependentTask) { /* Here we don't need a barrier (unlike the spreading), * since both tasks only construct vsites from particles, * or local vsites, not from non-local vsites. */ - construct_vsites_thread(x, dt, v, ip, tData.idTask.ilist, pbc_null); + construct_vsites_thread( + x, v, ip, tData.idTask.ilist, pbc_null); } } GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR } /* Now we can construct the vsites that might depend on other vsites */ - construct_vsites_thread(x, dt, v, ip, threadingInfo->threadDataNonLocalDependent().ilist, pbc_null); + construct_vsites_thread( + x, v, ip, threadingInfo->threadDataNonLocalDependent().ilist, pbc_null); } } -void VirtualSitesHandler::Impl::construct(ArrayRef x, real dt, ArrayRef v, const matrix box) const +void VirtualSitesHandler::Impl::construct(ArrayRef x, + ArrayRef v, + const matrix box, + VSiteOperation operation) const { - construct_vsites(&threadingInfo_, x, dt, v, iparams_, ilists_, domainInfo_, box); + switch (operation) + { + case VSiteOperation::Positions: + construct_vsites( + &threadingInfo_, x, v, iparams_, ilists_, domainInfo_, box); + break; + case VSiteOperation::Velocities: + construct_vsites( + &threadingInfo_, x, v, iparams_, ilists_, domainInfo_, box); + break; + case VSiteOperation::PositionsAndVelocities: + construct_vsites( + &threadingInfo_, x, v, iparams_, ilists_, domainInfo_, box); + break; + default: gmx_fatal(FARGS, "Unknown virtual site operation"); + } } -void VirtualSitesHandler::construct(ArrayRef x, real dt, ArrayRef v, const matrix box) const +void VirtualSitesHandler::construct(ArrayRef x, ArrayRef v, const matrix box, VSiteOperation operation) const { - impl_->construct(x, dt, v, box); + impl_->construct(x, v, box, operation); } void constructVirtualSites(ArrayRef x, ArrayRef ip, ArrayRef ilist) @@ -851,7 +1257,8 @@ void constructVirtualSites(ArrayRef x, ArrayRef ip, Array { // No PBC, no DD const DomainInfo domainInfo; - construct_vsites(nullptr, x, 0, {}, ip, ilist, domainInfo, nullptr); + construct_vsites( + nullptr, x, {}, ip, ilist, domainInfo, nullptr); } #ifndef DOXYGEN diff --git a/src/gromacs/mdlib/vsite.h b/src/gromacs/mdlib/vsite.h index 5fa0163d15..b8f694711d 100644 --- a/src/gromacs/mdlib/vsite.h +++ b/src/gromacs/mdlib/vsite.h @@ -80,6 +80,15 @@ static constexpr int c_ftypeVsiteEnd = F_VSITEN + 1; //! Type for storing PBC atom information for all vsite types in the system typedef std::array, c_ftypeVsiteEnd - c_ftypeVsiteStart> VsitePbc; +//! Whether we calculate vsite positions, velocities, or both +enum class VSiteOperation +{ + Positions, //!< Calculate only positions + Velocities, //!< Calculate only velocities + PositionsAndVelocities, //!< Calculate both positions and velocities + Count //!< The number of entries +}; + /*! \libinternal * \brief Class that handles construction of vsites and spreading of vsite forces */ @@ -99,12 +108,12 @@ public: /*! \brief Create positions of vsite atoms based for the local system * - * \param[in,out] x The coordinates - * \param[in] dt The time step - * \param[in,out] v When not empty, velocities for vsites are set as displacement/dt - * \param[in] box The box + * \param[in,out] x The coordinates + * \param[in,out] v The velocities, needed if operation requires it + * \param[in] box The box + * \param[in] operation Whether we calculate positions, velocities, or both */ - void construct(ArrayRef x, real dt, ArrayRef v, const matrix box) const; + void construct(ArrayRef x, ArrayRef v, const matrix box, VSiteOperation operation) const; //! Tells how to handle virial contributions due to virtual sites enum class VirialHandling : int diff --git a/src/gromacs/mdrun/md.cpp b/src/gromacs/mdrun/md.cpp index ab2d4b77b1..cf64bf19c6 100644 --- a/src/gromacs/mdrun/md.cpp +++ b/src/gromacs/mdrun/md.cpp @@ -551,7 +551,7 @@ void gmx::LegacySimulator::do_md() /* Set the velocities of vsites, shells and frozen atoms to zero */ for (i = 0; i < mdatoms->homenr; i++) { - if (mdatoms->ptype[i] == eptVSite || mdatoms->ptype[i] == eptShell) + if (mdatoms->ptype[i] == eptShell) { clear_rvec(v[i]); } @@ -949,11 +949,22 @@ void gmx::LegacySimulator::do_md() stateGpu->waitCoordinatesReadyOnHost(AtomLocality::Local); } + // We only need to calculate virtual velocities if we are writing them in the current step + const bool needVirtualVelocitiesThisStep = + (vsite != nullptr) + && (do_per_step(step, ir->nstvout) || checkpointHandler->isCheckpointingStep()); + if (vsite != nullptr) { // Virtual sites need to be updated before domain decomposition and forces are calculated wallcycle_start(wcycle, ewcVSITECONSTR); - vsite->construct(state->x, ir->delta_t, state->v, state->box); + // md-vv calculates virtual velocities once it has full-step real velocities + vsite->construct(state->x, + state->v, + state->box, + (!EI_VV(inputrec->eI) && needVirtualVelocitiesThisStep) + ? VSiteOperation::PositionsAndVelocities + : VSiteOperation::Positions); wallcycle_stop(wcycle, ewcVSITECONSTR); } @@ -1242,6 +1253,13 @@ void gmx::LegacySimulator::do_md() mdlog, fplog, wcycle); + if (vsite != nullptr && needVirtualVelocitiesThisStep) + { + // Positions were calculated earlier + wallcycle_start(wcycle, ewcVSITECONSTR); + vsite->construct(state->x, state->v, state->box, VSiteOperation::Velocities); + wallcycle_stop(wcycle, ewcVSITECONSTR); + } } /* ######## END FIRST UPDATE STEP ############## */ diff --git a/src/gromacs/mdrun/mimic.cpp b/src/gromacs/mdrun/mimic.cpp index 5e996d769b..d65d8b7898 100644 --- a/src/gromacs/mdrun/mimic.cpp +++ b/src/gromacs/mdrun/mimic.cpp @@ -472,7 +472,7 @@ void gmx::LegacySimulator::do_mimic() if (constructVsites) { wallcycle_start(wcycle, ewcVSITECONSTR); - vsite->construct(state->x, ir->delta_t, state->v, state->box); + vsite->construct(state->x, state->v, state->box, VSiteOperation::PositionsAndVelocities); wallcycle_stop(wcycle, ewcVSITECONSTR); } } diff --git a/src/gromacs/mdrun/minimize.cpp b/src/gromacs/mdrun/minimize.cpp index 07a9a4c5f6..a066a52655 100644 --- a/src/gromacs/mdrun/minimize.cpp +++ b/src/gromacs/mdrun/minimize.cpp @@ -911,7 +911,7 @@ void EnergyEvaluator::run(em_state_t* ems, rvec mu_tot, tensor vir, tensor pres, if (vsite) { - vsite->construct(ems->s.x, 1, {}, ems->s.box); + vsite->construct(ems->s.x, {}, ems->s.box, gmx::VSiteOperation::Positions); } if (DOMAINDECOMP(cr) && bNS) @@ -2022,7 +2022,7 @@ void LegacySimulator::do_lbfgs() if (vsite) { - vsite->construct(state_global->x, 1, {}, state_global->box); + vsite->construct(state_global->x, {}, state_global->box, VSiteOperation::Positions); } /* Call the force routine and some auxiliary (neighboursearching etc.) */ diff --git a/src/gromacs/mdrun/rerun.cpp b/src/gromacs/mdrun/rerun.cpp index 710e64c1bf..91857cb1c0 100644 --- a/src/gromacs/mdrun/rerun.cpp +++ b/src/gromacs/mdrun/rerun.cpp @@ -144,13 +144,11 @@ using gmx::VirtualSitesHandler; * \param[in,out] globalState The global state container * \param[in] constructVsites When true, vsite coordinates are constructed * \param[in] vsite Vsite setup, can be nullptr when \p constructVsites = false - * \param[in] timeStep Time step, used for constructing vsites */ static void prepareRerunState(const t_trxframe& rerunFrame, t_state* globalState, bool constructVsites, - const VirtualSitesHandler* vsite, - double timeStep) + const VirtualSitesHandler* vsite) { auto x = makeArrayRef(globalState->x); auto rerunX = arrayRefFromArray(reinterpret_cast(rerunFrame.x), globalState->natoms); @@ -161,7 +159,7 @@ static void prepareRerunState(const t_trxframe& rerunFrame, { GMX_ASSERT(vsite, "Need valid vsite for constructing vsites"); - vsite->construct(globalState->x, timeStep, globalState->v, globalState->box); + vsite->construct(globalState->x, globalState->v, globalState->box, gmx::VSiteOperation::PositionsAndVelocities); } } @@ -584,7 +582,7 @@ void gmx::LegacySimulator::do_rerun() "decomposition, " "use a single rank"); } - prepareRerunState(rerun_fr, state_global, constructVsites, vsite, ir->delta_t); + prepareRerunState(rerun_fr, state_global, constructVsites, vsite); } isLastStep = isLastStep || stopHandler->stoppingAfterCurrentStep(bNS); diff --git a/src/gromacs/mdrun/shellfc.cpp b/src/gromacs/mdrun/shellfc.cpp index 37522931b0..7521606b7b 100644 --- a/src/gromacs/mdrun/shellfc.cpp +++ b/src/gromacs/mdrun/shellfc.cpp @@ -1186,7 +1186,7 @@ void relax_shell_flexcon(FILE* fplog, { if (vsite) { - vsite->construct(pos[Min], inputrec->delta_t, v, box); + vsite->construct(pos[Min], v, box, gmx::VSiteOperation::PositionsAndVelocities); } if (nflexcon) diff --git a/src/gromacs/modularsimulator/statepropagatordata.cpp b/src/gromacs/modularsimulator/statepropagatordata.cpp index 866136a809..ce7a3c8618 100644 --- a/src/gromacs/modularsimulator/statepropagatordata.cpp +++ b/src/gromacs/modularsimulator/statepropagatordata.cpp @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2019,2020, by the GROMACS development team, led by + * Copyright (c) 2019,2020,2021, 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. @@ -148,7 +148,7 @@ StatePropagatorData::StatePropagatorData(int numAtoms, // Set the velocities of vsites, shells and frozen atoms to zero for (int i = 0; i < mdatoms->homenr; i++) { - if (mdatoms->ptype[i] == eptVSite || mdatoms->ptype[i] == eptShell) + if (mdatoms->ptype[i] == eptShell) { clear_rvec(v[i]); } diff --git a/src/programs/mdrun/tests/virtualsites.cpp b/src/programs/mdrun/tests/virtualsites.cpp index 9856a6937f..00f996beaf 100644 --- a/src/programs/mdrun/tests/virtualsites.cpp +++ b/src/programs/mdrun/tests/virtualsites.cpp @@ -134,8 +134,7 @@ public: velocities.emplace_back(frame.v().at(vSite.atomIdx)); } EXPECT_THAT(refPositions, Pointwise(RVecEq(tolerance), positions)); - // To be re-enabled in !979 once the velocity-handling is fixed - // EXPECT_THAT(refVelocities, Pointwise(RVecEq(tolerance), velocities)); + EXPECT_THAT(refVelocities, Pointwise(RVecEq(tolerance), velocities)); } } -- 2.22.0