#include "gromacs/utility/qsort_threadsafe.h"
#include "gromacs/pulling/pull_rotation.h"
#include "gromacs/mdlib/groupcoord.h"
+#include "gromacs/math/utilities.h"
static char *RotStr = {"Enforced rotation:"};
* However, since the atom is located directly on the pivot, this
* slab's contribution to the force on that atom will be zero
* anyway. Therefore, we directly move on to the next slab. */
- if (0 == norm(tmpvec2) )
+ if (gmx_numzero(norm(tmpvec2))) /* 0 == norm(xj - xcn) */
{
continue;
}
rvec_sub(yj0, ycn, yj0_ycn); /* yj0_ycn = yj0 - ycn */
- /* Rotate: */
- mvmul(erg->rotmat, yj0_ycn, tmpvec2); /* tmpvec2= Omega.(yj0-ycn) */
-
- /* Subtract the slab center from xj */
- rvec_sub(xj, xcn, xj_xcn); /* xj_xcn = xj - xcn */
-
- /* In rare cases, when an atom position coincides with a slab center
- * (xj_xcn == 0) we cannot compute the vector product for qjn.
+ /* In rare cases, when an atom position coincides with a reference slab
+ * center (yj0_ycn == 0) we cannot compute the normal vector qjn.
* However, since the atom is located directly on the pivot, this
* slab's contribution to the force on that atom will be zero
* anyway. Therefore, we directly move on to the next slab. */
- if (0 == norm(xj_xcn) )
+ if (gmx_numzero(norm(yj0_ycn))) /* 0 == norm(yj0 - ycn) */
{
continue;
}
+ /* Rotate: */
+ mvmul(erg->rotmat, yj0_ycn, tmpvec2); /* tmpvec2= Omega.(yj0-ycn) */
+
+ /* Subtract the slab center from xj */
+ rvec_sub(xj, xcn, xj_xcn); /* xj_xcn = xj - xcn */
+
/* Calculate qjn */
cprod(rotg->vec, tmpvec2, tmpvec); /* tmpvec= v x Omega.(yj0-ycn) */