Fixes #1431 (rotation/flex-t regression test failing on BG/Q)
In do_flex_lowlevel() we checked (by mistake!) for xj-xcn being
zero, although we need to check for yj0-ycn being zero, since
we use yj0-ycn in a cross product in the following lines of code.
I now also replaced the direct check (0 == norm(...)) by checking
what gmx_numzero(norm(...)) returns. The latter replacement
was also applied in the do_flex2_lowlevel() routine. Note that there
the check for small xj-xcn was and is actually correct.
Change-Id: I972b6d67a81e30f297db286cd2224f66753a20aa
#include "pull_rotation.h"
#include "gmx_sort.h"
#include "copyrite.h"
#include "pull_rotation.h"
#include "gmx_sort.h"
#include "copyrite.h"
static char *RotStr = {"Enforced rotation:"};
static char *RotStr = {"Enforced rotation:"};
/* Subtract the slab center from xj */
rvec_sub(xj, xcn, tmpvec2); /* tmpvec2 = xj - xcn */
/* Subtract the slab center from xj */
rvec_sub(xj, xcn, tmpvec2); /* tmpvec2 = xj - xcn */
/* In rare cases, when an atom position coincides with a slab center
/* In rare cases, when an atom position coincides with a slab center
- * (tmpvec2 == 0) we cannot compute the vector product for sjn.
- * However, since the atom is located directly on the pivot, this
- * slab's contribution to the force on that atom will be zero
+ * (tmpvec2 == 0) we cannot compute the vector product for sjn.
+ * 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. */
* anyway. Therefore, we directly move on to the next slab. */
- if ( 0 == norm(tmpvec2) )
+ if (gmx_numzero(norm(tmpvec2))) /* 0 == norm(xj - xcn) */
rvec_sub(yj0, ycn, yj0_ycn); /* yj0_ycn = yj0 - ycn */
rvec_sub(yj0, ycn, yj0_ycn); /* yj0_ycn = yj0 - ycn */
+ /* 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 (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 */
/* 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.
- * 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) )
- {
- continue;
- }
-
/* Calculate qjn */
cprod(rotg->vec, tmpvec2, tmpvec); /* tmpvec= v x Omega.(yj0-ycn) */
/* Calculate qjn */
cprod(rotg->vec, tmpvec2, tmpvec); /* tmpvec= v x Omega.(yj0-ycn) */