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:
.. math:: \mathbf{F}_i ~=~ \mathbf{F}_{s}
:label: eqnvsite1force
-- As a linear combination of two atoms
- (:numref:`Fig. %s <fig-vsites>` 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 <fig-vsites>` 2)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+- The weights are calculated as
.. math:: w_i = 1 - a ~,~~ w_j = a
:label: eqnvsitelin2atom
- 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 <fig-vsites>` 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 <fig-vsites>` 2fd)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+- The position is calculated as:
.. math:: \mathbf{r}_s ~=~ \mathbf{r}_i + a \frac{ \mathbf{r}_{ij} }
{ | \mathbf{r}_{ij} | }
\end{array}
:label: eqnvsite2fdforce
-- As a linear combination of three atoms
- (:numref:`Fig. %s <fig-vsites>` 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 <fig-vsites>` 3)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+- The weights are calculated as:
.. math:: w_i = 1 - a - b ~,~~ w_j = a ~,~~ w_k = b
:label: eqnvsitelin3atom
- 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 <fig-vsites>` 3fd):
+In the plane of three atoms, with a fixed distance (:numref:`Fig. %s <fig-vsites>` 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
\end{array}
:label: eqnvsiteplane3atomforce
-- In the plane of three atoms, with a fixed angle and
- distance (:numref:`Fig. %s <fig-vsites>` 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 <fig-vsites>` 3fad)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+- The position is calculated as:
.. math:: \mathbf{r}_s ~=~ \mathbf{r}_i +
d \cos \theta \frac{\mathbf{r}_{ij}}{ | \mathbf{r}_{ij} | } +
\end{array}
:label: eqnvsite2fadFforce
-- As a non-linear combination of three atoms, out of
- plane (:numref:`Fig. %s <fig-vsites>` 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 <fig-vsites>` 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})
\end{array}
:label: eqnvsitenonlin3atomforce
-- From four atoms, with a fixed distance, see
- separate :numref:`Fig. %s <fig-vsite4fdn>`. 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 <fig-vsite4fdn>`
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+- 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:
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}`
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
* 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.
/*! \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<RVec> x, real dt, ArrayRef<RVec> v, const matrix box) const;
+ void construct(ArrayRef<RVec> x, ArrayRef<RVec> v, const matrix box, VSiteOperation operation) const;
/*! \brief Spread the force operating on the vsite atoms on the surrounding atoms.
*
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<VSiteCalculatePosition calculatePosition, VSiteCalculateVelocity calculateVelocity>
+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<VSiteCalculatePosition calculatePosition, VSiteCalculateVelocity calculateVelocity>
+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<VSiteCalculatePosition calculatePosition, VSiteCalculateVelocity calculateVelocity>
+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<VSiteCalculatePosition calculatePosition, VSiteCalculateVelocity calculateVelocity>
+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<VSiteCalculatePosition calculatePosition, VSiteCalculateVelocity calculateVelocity>
+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);
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<VSiteCalculatePosition calculatePosition, VSiteCalculateVelocity calculateVelocity>
+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<VSiteCalculatePosition calculatePosition, VSiteCalculateVelocity calculateVelocity>
+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;
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<VSiteCalculatePosition calculatePosition, VSiteCalculateVelocity calculateVelocity>
static void constr_vsite4FD(const rvec xi,
const rvec xj,
const rvec xk,
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;
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<VSiteCalculatePosition calculatePosition, VSiteCalculateVelocity calculateVelocity>
static void constr_vsite4FDN(const rvec xi,
const rvec xj,
const rvec xk,
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;
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<const t_iparams> ip, ArrayRef<RVec> x, const t_pbc* pbc)
+template<VSiteCalculatePosition calculatePosition, VSiteCalculateVelocity calculateVelocity>
+static int constr_vsiten(const t_iatom* ia,
+ ArrayRef<const t_iparams> ip,
+ ArrayRef<RVec> x,
+ const t_pbc* pbc,
+ ArrayRef<RVec> 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
/*! \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<VSiteCalculatePosition calculatePosition, VSiteCalculateVelocity calculateVelocity>
static void construct_vsites_thread(ArrayRef<RVec> x,
- const real dt,
ArrayRef<RVec> v,
ArrayRef<const t_iparams> ip,
ArrayRef<const InteractionList> 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 */
real b1, c1;
switch (ftype)
{
- case F_VSITE1: constr_vsite1(x[ai], x[avsite]); break;
+ case F_VSITE1:
+ constr_vsite1<calculatePosition, calculateVelocity>(
+ 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<calculatePosition, calculateVelocity>(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<calculatePosition, calculateVelocity>(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<calculatePosition, calculateVelocity>(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<calculatePosition, calculateVelocity>(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<calculatePosition, calculateVelocity>(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<calculatePosition, calculateVelocity>(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];
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<calculatePosition, calculateVelocity>(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];
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<calculatePosition, calculateVelocity>(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<calculatePosition, calculateVelocity>(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__);
}
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;
*
* \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<VSiteCalculatePosition calculatePosition, VSiteCalculateVelocity calculateVelocity>
static void construct_vsites(const ThreadingInfo* threadingInfo,
ArrayRef<RVec> x,
- real dt,
ArrayRef<RVec> v,
ArrayRef<const t_iparams> ip,
ArrayRef<const InteractionList> ilist,
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<calculatePosition, calculateVelocity>(x, v, ip, ilist, pbc_null);
}
else
{
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<calculatePosition, calculateVelocity>(
+ 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<calculatePosition, calculateVelocity>(
+ 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<calculatePosition, calculateVelocity>(
+ x, v, ip, threadingInfo->threadDataNonLocalDependent().ilist, pbc_null);
}
}
-void VirtualSitesHandler::Impl::construct(ArrayRef<RVec> x, real dt, ArrayRef<RVec> v, const matrix box) const
+void VirtualSitesHandler::Impl::construct(ArrayRef<RVec> x,
+ ArrayRef<RVec> v,
+ const matrix box,
+ VSiteOperation operation) const
{
- construct_vsites(&threadingInfo_, x, dt, v, iparams_, ilists_, domainInfo_, box);
+ switch (operation)
+ {
+ case VSiteOperation::Positions:
+ construct_vsites<VSiteCalculatePosition::Yes, VSiteCalculateVelocity::No>(
+ &threadingInfo_, x, v, iparams_, ilists_, domainInfo_, box);
+ break;
+ case VSiteOperation::Velocities:
+ construct_vsites<VSiteCalculatePosition::No, VSiteCalculateVelocity::Yes>(
+ &threadingInfo_, x, v, iparams_, ilists_, domainInfo_, box);
+ break;
+ case VSiteOperation::PositionsAndVelocities:
+ construct_vsites<VSiteCalculatePosition::Yes, VSiteCalculateVelocity::Yes>(
+ &threadingInfo_, x, v, iparams_, ilists_, domainInfo_, box);
+ break;
+ default: gmx_fatal(FARGS, "Unknown virtual site operation");
+ }
}
-void VirtualSitesHandler::construct(ArrayRef<RVec> x, real dt, ArrayRef<RVec> v, const matrix box) const
+void VirtualSitesHandler::construct(ArrayRef<RVec> x, ArrayRef<RVec> v, const matrix box, VSiteOperation operation) const
{
- impl_->construct(x, dt, v, box);
+ impl_->construct(x, v, box, operation);
}
void constructVirtualSites(ArrayRef<RVec> x, ArrayRef<const t_iparams> ip, ArrayRef<const InteractionList> ilist)
{
// No PBC, no DD
const DomainInfo domainInfo;
- construct_vsites(nullptr, x, 0, {}, ip, ilist, domainInfo, nullptr);
+ construct_vsites<VSiteCalculatePosition::Yes, VSiteCalculateVelocity::No>(
+ nullptr, x, {}, ip, ilist, domainInfo, nullptr);
}
#ifndef DOXYGEN
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff