From e16a45d7cf4668ec98b07b1db00bf5da0d973dc3 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).

Backported from master with minor changes, fixes #1611.  Kept commit
message the same; the second point will be only relevant for master.

Change-Id: Ib7602fa49a1b8f7882a63843322786b3e51e8e32
(cherry-picked from b3e2e82 in master)
---
 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 24a66ddbda..599d421881 100644
--- a/src/gromacs/selection/nbsearch.cpp
+++ b/src/gromacs/selection/nbsearch.cpp
@@ -131,10 +131,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.
          *
@@ -423,25 +423,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
@@ -494,15 +521,8 @@ void AnalysisNeighborhoodSearchImpl::init(
 
         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);
         }
     }
@@ -543,12 +563,13 @@ void AnalysisNeighborhoodPairSearchImpl::reset(int testIndex)
     testIndex_ = testIndex;
     if (testIndex_ >= 0 && testIndex_ < static_cast<int>(testPositions_.size()))
     {
-        copy_rvec(testPositions_[testIndex_], xtest_);
         if (search_.bGrid_)
         {
-            put_atoms_in_triclinic_unitcell(ecenterTRIC, search_.pbc_->box,
-                                            1, &xtest_);
-            search_.mapPointToGridCell(xtest_, testcell_);
+            search_.mapPointToGridCell(testPositions_[testIndex], testcell_, xtest_);
+        }
+        else
+        {
+            copy_rvec(testPositions_[testIndex_], xtest_);
         }
     }
     previ_     = -1;
-- 
2.11.4.GIT