From: Carsten Kutzner <ckutzne@gwdg.de>
Date: Thu, 23 May 2013 16:23:41 +0000 (+0200)
Subject: Fix for a rare condition in flexible rotation potentials.
X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=f762a4f8439aadb863ea54ddc342871d2122057c;p=alexxy%2Fgromacs.git

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 8fbe2272fa..e51d4d0871 100644
--- 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) */