+ 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);
+ }