From b207685897ed4ab3cea11008c29d4b6005cf324d Mon Sep 17 00:00:00 2001 From: Carsten Kutzner Date: Wed, 12 Mar 2014 16:14:04 +0100 Subject: [PATCH] Avoid cross product with zero vector in rotational pulling. 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 --- src/mdlib/pull_rotation.c | 31 ++++++++++++++++--------------- 1 file changed, 16 insertions(+), 15 deletions(-) diff --git a/src/mdlib/pull_rotation.c b/src/mdlib/pull_rotation.c index 1039daaeb1..54f558d0bf 100644 --- a/src/mdlib/pull_rotation.c +++ b/src/mdlib/pull_rotation.c @@ -64,6 +64,7 @@ #include "pull_rotation.h" #include "gmx_sort.h" #include "copyrite.h" +#include "maths.h" static char *RotStr = {"Enforced rotation:"}; @@ -1967,13 +1968,13 @@ static real do_flex2_lowlevel( /* 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 - * (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. */ - if ( 0 == norm(tmpvec2) ) + if (gmx_numzero(norm(tmpvec2))) /* 0 == norm(xj - xcn) */ { continue; } @@ -2201,22 +2202,22 @@ static real do_flex_lowlevel( 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 */ - /* 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) */ -- 2.22.0