From c929b623469f26e9519059e2cde445f9bb3cfc49 Mon Sep 17 00:00:00 2001
From: Berk Hess <hess@kth.se>
Date: Thu, 31 Jul 2014 13:08:34 +0200
Subject: [PATCH] Avoid division by 0 in nbnxn_search

Note that this issue did not affect any results.

Change-Id: I3c454783484c416b50808b9512e14072cb749784
---
 src/gromacs/mdlib/nbnxn_search.c | 20 +++++++++++++++++++-
 src/gromacs/mdlib/nbnxn_search.h |  2 +-
 2 files changed, 20 insertions(+), 2 deletions(-)

diff --git a/src/gromacs/mdlib/nbnxn_search.c b/src/gromacs/mdlib/nbnxn_search.c
index 1bc3f05847..dc0b02013f 100644
--- a/src/gromacs/mdlib/nbnxn_search.c
+++ b/src/gromacs/mdlib/nbnxn_search.c
@@ -419,6 +419,12 @@ static real grid_atom_density(int n, rvec corner0, rvec corner1)
 {
     rvec size;
 
+    if (n == 0)
+    {
+        /* To avoid zero density we use a minimum of 1 atom */
+        n = 1;
+    }
+
     rvec_sub(corner1, corner0, size);
 
     return n/(size[XX]*size[YY]*size[ZZ]);
@@ -439,6 +445,8 @@ static int set_grid_size_xy(const nbnxn_search_t nbs,
 
     if (n > grid->na_sc)
     {
+        assert(atom_density > 0);
+
         /* target cell length */
         if (grid->bSimple)
         {
@@ -1814,7 +1822,8 @@ void nbnxn_put_on_grid(nbnxn_search_t nbs,
         nbs->ePBC = ePBC;
         copy_mat(box, nbs->box);
 
-        if (atom_density >= 0)
+        /* Avoid zero density */
+        if (atom_density > 0)
         {
             grid->atom_density = atom_density;
         }
@@ -1830,12 +1839,21 @@ void nbnxn_put_on_grid(nbnxn_search_t nbs,
          * for the local atoms (dd_zone=0).
          */
         nbs->natoms_nonlocal = a1 - nmoved;
+
+        if (debug)
+        {
+            fprintf(debug, "natoms_local = %5d atom_density = %5.1f\n",
+                    nbs->natoms_local, grid->atom_density);
+        }
     }
     else
     {
         nbs->natoms_nonlocal = max(nbs->natoms_nonlocal, a1);
     }
 
+    /* We always use the home zone (grid[0]) for setting the cell size,
+     * since determining densities for non-local zones is difficult.
+     */
     nc_max_grid = set_grid_size_xy(nbs, grid,
                                    dd_zone, n-nmoved, corner0, corner1,
                                    nbs->grid[0].atom_density);
diff --git a/src/gromacs/mdlib/nbnxn_search.h b/src/gromacs/mdlib/nbnxn_search.h
index 12576a096c..6b3ab7c8d2 100644
--- a/src/gromacs/mdlib/nbnxn_search.h
+++ b/src/gromacs/mdlib/nbnxn_search.h
@@ -65,7 +65,7 @@ void nbnxn_init_search(nbnxn_search_t    * nbs_ptr,
 /* Put the atoms on the pair search grid.
  * Only atoms a0 to a1 in x are put on the grid.
  * The atom_density is used to determine the grid size.
- * When atom_density=-1, the density is determined from a1-a0 and the corners.
+ * When atom_density<=0, the density is determined from a1-a0 and the corners.
  * With domain decomposition part of the n particles might have migrated,
  * but have not been removed yet. This count is given by nmoved.
  * When move[i] < 0 particle i has migrated and will not be put on the grid.
-- 
2.22.0