Fix genion -nq and -pq when using -neutral

Currently genion can't neutralize (-neutral) systems when -nq/-pq
differs from -1/+1. This commit fixes the problem.

Change-Id: If2a82ace079555d8a684ea876daa947677d5e442

diff --git a/src/tools/gmx_genion.c b/src/tools/gmx_genion.c
index 95023713e4a0933bd68ae478468eae9f3406fb01..7c05a401b390ee83e5bf143db3a7827feb7f1080 100644 (file)
--- a/src/tools/gmx_genion.c
+++ b/src/tools/gmx_genion.c
@@ -63,6 +63,18 @@
 #include "mtop_util.h"
 #include "gmx_ana.h"
 
+static int greatest_common_divisor(int p, int q)
+{
+    int tmp;
+    while (q != 0)
+    {
+        tmp = q;
+        q = p % q;
+        p = tmp;
+    }
+    return p;
+}
+
 static void insert_ion(int nsa, int *nwater,
                        gmx_bool bSet[], int repl[], atom_id index[],
                        real pot[], rvec x[], t_pbc *pbc,
@@ -509,17 +521,28 @@ int gmx_genion(int argc, char *argv[])
     if (bNeutral)
     {
         int qdelta = p_num*p_q + n_num*n_q + iqtot;
+
+        /* Check if the system is neutralizable
+         * is (qdelta == p_q*p_num + n_q*n_num) solvable for p_num and n_num? */
+        int gcd = greatest_common_divisor(n_q, p_q);
+        if ((qdelta % gcd) != 0)
+        {
+            gmx_fatal(FARGS, "Can't neutralize this system using -nq %d and"
+                    " -pq %d.\n", n_q, p_q);
+        }
+        
         while (qdelta != 0)
         {
-            if (qdelta < 0)
+            while (qdelta < 0)
             {
-                p_num  += abs(qdelta/p_q);
+                p_num++;
+                qdelta += p_q;
             }
-            else if (qdelta > 0)
+            while (qdelta > 0)
             {
-                n_num  += abs(qdelta/n_q);
+                n_num++;
+                qdelta += n_q;
             }
-            qdelta  = p_num*p_q + n_num*n_q + iqtot;
         }
     }