Fix for a rare condition in flexible rotation potentials.

Change-Id: Icc6c7fceaba3fa53f3f2d627c3c205469fc9d710

diff --git a/src/mdlib/pull_rotation.c b/src/mdlib/pull_rotation.c
index 8fbe2272fa094b7e92ba47d48427406253cbdbbf..e51d4d0871e018197ea9fb0492aaa38b94fa3edc 100644 (file)
--- a/src/mdlib/pull_rotation.c
+++ b/src/mdlib/pull_rotation.c
@@ -1967,6 +1967,16 @@ 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 
+             * anyway. Therefore, we directly move on to the next slab.       */
+            if ( 0 == norm(tmpvec2) )
+            {
+                continue;
+            }
 
             /* Calculate sjn */
             cprod(rotg->vec, tmpvec2, tmpvec);    /* tmpvec = v x (xj - xcn)  */
@@ -2197,6 +2207,16 @@ static real do_flex_lowlevel(
             /* 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) */