From b3e2e828fbd7490662f2734834848d33e7fd6c46 Mon Sep 17 00:00:00 2001
From: Teemu Murtola <teemu.murtola@gmail.com>
Date: Sat, 23 Aug 2014 06:27:24 +0300
Subject: [PATCH] Improve analysis nbsearch grid mapping

Now the analysis neighborhood search implements its own version of
put_atoms_in_triclinic_unitcell().  While computing the index of the
correct grid cell, it is relatively easy to produce also the coordinates
that lay within that cell instead of using a separate call.  This
provides two benefits:
 - It avoids rare rounding problems if put_atoms_in_triclinic_unitcell()
   would put the atom right at the edge of the box, but the mapping code
   would consider it outside the box, causing out-of-range grid cell
   index to be generated.
 - It allows to customize the grid mapping more freely (e.g., to create
   grids that are not periodic).

Change-Id: Ib7602fa49a1b8f7882a63843322786b3e51e8e32
---
 src/gromacs/selection/nbsearch.cpp | 61 ++++++++++++++++++++----------
 1 file changed, 41 insertions(+), 20 deletions(-)

diff --git a/src/gromacs/selection/nbsearch.cpp b/src/gromacs/selection/nbsearch.cpp
index 758d97aa12..0b50214a4d 100644
--- a/src/gromacs/selection/nbsearch.cpp
+++ b/src/gromacs/selection/nbsearch.cpp
@@ -141,10 +141,10 @@ class AnalysisNeighborhoodSearchImpl
          *
          * \param[in]  x    Point to map.
          * \param[out] cell Indices of the grid cell in which \p x lies.
-         *
-         * \p x should be within the triclinic unit cell.
+         * \param[out] xout Coordinates to use
+         *     (will be within the triclinic unit cell).
          */
-        void mapPointToGridCell(const rvec x, ivec cell) const;
+        void mapPointToGridCell(const rvec x, ivec cell, rvec xout) const;
         /*! \brief
          * Calculates linear index of a grid cell.
          *
@@ -441,25 +441,52 @@ bool AnalysisNeighborhoodSearchImpl::initGrid(const t_pbc &pbc)
 }
 
 void AnalysisNeighborhoodSearchImpl::mapPointToGridCell(const rvec x,
-                                                        ivec       cell) const
+                                                        ivec       cell,
+                                                        rvec       xout) const
 {
+    rvec xtmp;
+    copy_rvec(x, xtmp);
     if (bTric_)
     {
         rvec tx;
-
-        tmvmul_ur0(recipcell_, x, tx);
+        tmvmul_ur0(recipcell_, xtmp, tx);
         for (int dd = 0; dd < DIM; ++dd)
         {
-            cell[dd] = static_cast<int>(tx[dd]);
+            const int cellCount = ncelldim_[dd];
+            int       cellIndex = static_cast<int>(floor(tx[dd]));
+            while (cellIndex < 0)
+            {
+                cellIndex += cellCount;
+                rvec_add(xtmp, pbc_.box[dd], xtmp);
+            }
+            while (cellIndex >= cellCount)
+            {
+                cellIndex -= cellCount;
+                rvec_sub(xtmp, pbc_.box[dd], xtmp);
+            }
+            cell[dd] = cellIndex;
         }
     }
     else
     {
         for (int dd = 0; dd < DIM; ++dd)
         {
-            cell[dd] = static_cast<int>(x[dd] * recipcell_[dd][dd]);
+            const int cellCount = ncelldim_[dd];
+            int       cellIndex = static_cast<int>(floor(xtmp[dd] * recipcell_[dd][dd]));
+            while (cellIndex < 0)
+            {
+                cellIndex += cellCount;
+                xtmp[dd]  += pbc_.box[dd][dd];
+            }
+            while (cellIndex >= cellCount)
+            {
+                cellIndex -= cellCount;
+                xtmp[dd]  -= pbc_.box[dd][dd];
+            }
+            cell[dd] = cellIndex;
         }
     }
+    copy_rvec(xtmp, xout);
 }
 
 int AnalysisNeighborhoodSearchImpl::getGridCellIndex(const ivec cell) const
@@ -536,17 +563,10 @@ void AnalysisNeighborhoodSearchImpl::init(
         }
         xref_ = xref_alloc_;
 
-        for (int i = 0; i < nref_; ++i)
-        {
-            copy_rvec(positions.x_[i], xref_alloc_[i]);
-        }
-        put_atoms_in_triclinic_unitcell(ecenterTRIC, pbc_.box,
-                                        nref_, xref_alloc_);
         for (int i = 0; i < nref_; ++i)
         {
             ivec refcell;
-
-            mapPointToGridCell(xref_[i], refcell);
+            mapPointToGridCell(positions.x_[i], refcell, xref_alloc_[i]);
             addToGridCell(refcell, i);
         }
     }
@@ -578,10 +598,11 @@ void AnalysisNeighborhoodPairSearchImpl::reset(int testIndex)
         copy_rvec(testPositions_[testIndex_], xtest_);
         if (search_.bGrid_)
         {
-            put_atoms_in_triclinic_unitcell(ecenterTRIC,
-                                            const_cast<rvec *>(search_.pbc_.box),
-                                            1, &xtest_);
-            search_.mapPointToGridCell(xtest_, testcell_);
+            search_.mapPointToGridCell(testPositions_[testIndex], testcell_, xtest_);
+        }
+        else
+        {
+            copy_rvec(testPositions_[testIndex_], xtest_);
         }
         if (search_.excls_ != NULL)
         {
-- 
2.22.0