src/gromacs/selection/tests/nbsearch.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2013,2014, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \internal \file
  36  * \brief
  37  * Tests selection neighborhood searching.
  38  *
  39  * \todo
  40  * Increase coverage of these tests for different corner cases: other PBC cases
  41  * than full 3D, large cutoffs (larger than half the box size), etc.
  42  * At least some of these probably don't work correctly.
  43  *
  44  * \author Teemu Murtola <teemu.murtola@gmail.com>
  45  * \ingroup module_selection
  46  */
  47 #include "gromacs/selection/nbsearch.h"
  48
  49 #include <gtest/gtest.h>
  50
  51 #include <cmath>
  52
  53 #include <algorithm>
  54 #include <limits>
  55 #include <vector>
  56
  57 #include "gromacs/math/vec.h"
  58 #include "gromacs/pbcutil/pbc.h"
  59 #include "gromacs/random/random.h"
  60 #include "gromacs/utility/smalloc.h"
  61
  62 #include "testutils/testasserts.h"
  63
  64 namespace
  65 {
  66
  67 /********************************************************************
  68  * NeighborhoodSearchTestData
  69  */
  70
  71 class NeighborhoodSearchTestData
  72 {
  73     public:
  74         struct RefPair
  75         {
  76             RefPair(int refIndex, real distance)
  77                 : refIndex(refIndex), distance(distance), bFound(false)
  78             {
  79             }
  80
  81             bool operator<(const RefPair &other) const
  82             {
  83                 return refIndex < other.refIndex;
  84             }
  85
  86             int                 refIndex;
  87             real                distance;
  88             bool                bFound;
  89         };
  90
  91         struct TestPosition
  92         {
  93             TestPosition() : refMinDist(0.0), refNearestPoint(-1)
  94             {
  95                 clear_rvec(x);
  96             }
  97             explicit TestPosition(const rvec x)
  98                 : refMinDist(0.0), refNearestPoint(-1)
  99             {
 100                 copy_rvec(x, this->x);
 101             }
 102
 103             rvec                 x;
 104             real                 refMinDist;
 105             int                  refNearestPoint;
 106             std::vector<RefPair> refPairs;
 107         };
 108
 109         typedef std::vector<TestPosition> TestPositionList;
 110
 111         NeighborhoodSearchTestData(int seed, real cutoff);
 112         ~NeighborhoodSearchTestData();
 113
 114         gmx::AnalysisNeighborhoodPositions refPositions() const
 115         {
 116             return gmx::AnalysisNeighborhoodPositions(refPos_, refPosCount_);
 117         }
 118         gmx::AnalysisNeighborhoodPositions testPositions() const
 119         {
 120             if (testPos_ == NULL)
 121             {
 122                 snew(testPos_, testPositions_.size());
 123                 for (size_t i = 0; i < testPositions_.size(); ++i)
 124                 {
 125                     copy_rvec(testPositions_[i].x, testPos_[i]);
 126                 }
 127             }
 128             return gmx::AnalysisNeighborhoodPositions(testPos_,
 129                                                       testPositions_.size());
 130         }
 131         gmx::AnalysisNeighborhoodPositions testPosition(int index) const
 132         {
 133             return testPositions().selectSingleFromArray(index);
 134         }
 135
 136         void addTestPosition(const rvec x)
 137         {
 138             GMX_RELEASE_ASSERT(testPos_ == NULL,
 139                                "Cannot add positions after testPositions() call");
 140             testPositions_.push_back(TestPosition(x));
 141         }
 142         void generateRandomPosition(rvec x);
 143         void generateRandomRefPositions(int count);
 144         void generateRandomTestPositions(int count);
 145         void computeReferences(t_pbc *pbc)
 146         {
 147             computeReferencesInternal(pbc, false);
 148         }
 149         void computeReferencesXY(t_pbc *pbc)
 150         {
 151             computeReferencesInternal(pbc, true);
 152         }
 153
 154         bool containsPair(int testIndex, const RefPair &pair) const
 155         {
 156             const std::vector<RefPair>          &refPairs = testPositions_[testIndex].refPairs;
 157             std::vector<RefPair>::const_iterator foundRefPair
 158                 = std::lower_bound(refPairs.begin(), refPairs.end(), pair);
 159             if (foundRefPair == refPairs.end() || foundRefPair->refIndex != pair.refIndex)
 160             {
 161                 return false;
 162             }
 163             return true;
 164         }
 165
 166         gmx_rng_t                        rng_;
 167         real                             cutoff_;
 168         matrix                           box_;
 169         t_pbc                            pbc_;
 170         int                              refPosCount_;
 171         rvec                            *refPos_;
 172         TestPositionList                 testPositions_;
 173
 174     private:
 175         void computeReferencesInternal(t_pbc *pbc, bool bXY);
 176
 177         mutable rvec                    *testPos_;
 178 };
 179
 180 //! Shorthand for a collection of reference pairs.
 181 typedef std::vector<NeighborhoodSearchTestData::RefPair> RefPairList;
 182
 183 NeighborhoodSearchTestData::NeighborhoodSearchTestData(int seed, real cutoff)
 184     : rng_(NULL), cutoff_(cutoff), refPosCount_(0), refPos_(NULL), testPos_(NULL)
 185 {
 186     // TODO: Handle errors.
 187     rng_ = gmx_rng_init(seed);
 188     clear_mat(box_);
 189     set_pbc(&pbc_, epbcNONE, box_);
 190 }
 191
 192 NeighborhoodSearchTestData::~NeighborhoodSearchTestData()
 193 {
 194     if (rng_ != NULL)
 195     {
 196         gmx_rng_destroy(rng_);
 197     }
 198     sfree(refPos_);
 199     sfree(testPos_);
 200 }
 201
 202 void NeighborhoodSearchTestData::generateRandomPosition(rvec x)
 203 {
 204     rvec fx;
 205     fx[XX] = gmx_rng_uniform_real(rng_);
 206     fx[YY] = gmx_rng_uniform_real(rng_);
 207     fx[ZZ] = gmx_rng_uniform_real(rng_);
 208     mvmul(box_, fx, x);
 209     // Add a small displacement to allow positions outside the box
 210     x[XX] += 0.2 * gmx_rng_uniform_real(rng_) - 0.1;
 211     x[YY] += 0.2 * gmx_rng_uniform_real(rng_) - 0.1;
 212     x[ZZ] += 0.2 * gmx_rng_uniform_real(rng_) - 0.1;
 213 }
 214
 215 void NeighborhoodSearchTestData::generateRandomRefPositions(int count)
 216 {
 217     refPosCount_ = count;
 218     snew(refPos_, refPosCount_);
 219     for (int i = 0; i < refPosCount_; ++i)
 220     {
 221         generateRandomPosition(refPos_[i]);
 222     }
 223 }
 224
 225 void NeighborhoodSearchTestData::generateRandomTestPositions(int count)
 226 {
 227     testPositions_.reserve(count);
 228     for (int i = 0; i < count; ++i)
 229     {
 230         rvec x;
 231         generateRandomPosition(x);
 232         addTestPosition(x);
 233     }
 234 }
 235
 236 void NeighborhoodSearchTestData::computeReferencesInternal(t_pbc *pbc, bool bXY)
 237 {
 238     real cutoff = cutoff_;
 239     if (cutoff <= 0)
 240     {
 241         cutoff = std::numeric_limits<real>::max();
 242     }
 243     TestPositionList::iterator i;
 244     for (i = testPositions_.begin(); i != testPositions_.end(); ++i)
 245     {
 246         i->refMinDist      = cutoff;
 247         i->refNearestPoint = -1;
 248         i->refPairs.clear();
 249         for (int j = 0; j < refPosCount_; ++j)
 250         {
 251             rvec dx;
 252             if (pbc != NULL)
 253             {
 254                 pbc_dx(pbc, i->x, refPos_[j], dx);
 255             }
 256             else
 257             {
 258                 rvec_sub(i->x, refPos_[j], dx);
 259             }
 260             // TODO: This may not work intuitively for 2D with the third box
 261             // vector not parallel to the Z axis, but neither does the actual
 262             // neighborhood search.
 263             const real dist =
 264                 !bXY ? norm(dx) : sqrt(sqr(dx[XX]) + sqr(dx[YY]));
 265             if (dist < i->refMinDist)
 266             {
 267                 i->refMinDist      = dist;
 268                 i->refNearestPoint = j;
 269             }
 270             if (dist <= cutoff)
 271             {
 272                 RefPair pair(j, dist);
 273                 GMX_RELEASE_ASSERT(i->refPairs.empty() || i->refPairs.back() < pair,
 274                                    "Reference pairs should be generated in sorted order");
 275                 i->refPairs.push_back(pair);
 276             }
 277         }
 278     }
 279 }
 280
 281 /********************************************************************
 282  * NeighborhoodSearchTest
 283  */
 284
 285 class NeighborhoodSearchTest : public ::testing::Test
 286 {
 287     public:
 288         void testIsWithin(gmx::AnalysisNeighborhoodSearch  *search,
 289                           const NeighborhoodSearchTestData &data);
 290         void testMinimumDistance(gmx::AnalysisNeighborhoodSearch  *search,
 291                                  const NeighborhoodSearchTestData &data);
 292         void testNearestPoint(gmx::AnalysisNeighborhoodSearch  *search,
 293                               const NeighborhoodSearchTestData &data);
 294         void testPairSearch(gmx::AnalysisNeighborhoodSearch  *search,
 295                             const NeighborhoodSearchTestData &data);
 296
 297         gmx::AnalysisNeighborhood        nb_;
 298 };
 299
 300 void NeighborhoodSearchTest::testIsWithin(
 301         gmx::AnalysisNeighborhoodSearch  *search,
 302         const NeighborhoodSearchTestData &data)
 303 {
 304     NeighborhoodSearchTestData::TestPositionList::const_iterator i;
 305     for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
 306     {
 307         const bool bWithin = (i->refMinDist <= data.cutoff_);
 308         EXPECT_EQ(bWithin, search->isWithin(i->x))
 309         << "Distance is " << i->refMinDist;
 310     }
 311 }
 312
 313 void NeighborhoodSearchTest::testMinimumDistance(
 314         gmx::AnalysisNeighborhoodSearch  *search,
 315         const NeighborhoodSearchTestData &data)
 316 {
 317     NeighborhoodSearchTestData::TestPositionList::const_iterator i;
 318     for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
 319     {
 320         const real refDist = i->refMinDist;
 321         EXPECT_REAL_EQ_TOL(refDist, search->minimumDistance(i->x),
 322                            gmx::test::ulpTolerance(20));
 323     }
 324 }
 325
 326 void NeighborhoodSearchTest::testNearestPoint(
 327         gmx::AnalysisNeighborhoodSearch  *search,
 328         const NeighborhoodSearchTestData &data)
 329 {
 330     NeighborhoodSearchTestData::TestPositionList::const_iterator i;
 331     for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
 332     {
 333         const gmx::AnalysisNeighborhoodPair pair = search->nearestPoint(i->x);
 334         if (pair.isValid())
 335         {
 336             EXPECT_EQ(i->refNearestPoint, pair.refIndex());
 337             EXPECT_EQ(0, pair.testIndex());
 338             EXPECT_REAL_EQ_TOL(i->refMinDist, sqrt(pair.distance2()),
 339                                gmx::test::ulpTolerance(64));
 340         }
 341         else
 342         {
 343             EXPECT_EQ(i->refNearestPoint, -1);
 344         }
 345     }
 346 }
 347
 348 void NeighborhoodSearchTest::testPairSearch(
 349         gmx::AnalysisNeighborhoodSearch  *search,
 350         const NeighborhoodSearchTestData &data)
 351 {
 352     NeighborhoodSearchTestData::TestPositionList::const_iterator i;
 353     // TODO: Test also searching all the test positions in a single search;
 354     // currently the implementation just contains this loop, though, but in
 355     // the future that may trigger a different code path.
 356     for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
 357     {
 358         RefPairList                         refPairs = i->refPairs;
 359         gmx::AnalysisNeighborhoodPairSearch pairSearch
 360             = search->startPairSearch(i->x);
 361         gmx::AnalysisNeighborhoodPair       pair;
 362         while (pairSearch.findNextPair(&pair))
 363         {
 364             EXPECT_EQ(0, pair.testIndex());
 365             NeighborhoodSearchTestData::RefPair searchPair(pair.refIndex(),
 366                                                            sqrt(pair.distance2()));
 367             RefPairList::iterator               foundRefPair
 368                 = std::lower_bound(refPairs.begin(), refPairs.end(), searchPair);
 369             if (foundRefPair == refPairs.end() || foundRefPair->refIndex != pair.refIndex())
 370             {
 371                 ADD_FAILURE()
 372                 << "Expected: Position " << pair.refIndex()
 373                 << " is within cutoff.\n"
 374                 << "  Actual: It is not.";
 375             }
 376             else if (foundRefPair->bFound)
 377             {
 378                 ADD_FAILURE()
 379                 << "Expected: Position " << pair.refIndex()
 380                 << " is returned only once.\n"
 381                 << "  Actual: It is returned multiple times.";
 382             }
 383             else
 384             {
 385                 foundRefPair->bFound = true;
 386                 EXPECT_REAL_EQ_TOL(foundRefPair->distance, searchPair.distance,
 387                                    gmx::test::ulpTolerance(64));
 388             }
 389         }
 390         for (RefPairList::const_iterator j = refPairs.begin(); j != refPairs.end(); ++j)
 391         {
 392             if (!j->bFound)
 393             {
 394                 ADD_FAILURE()
 395                 << "Expected: All pairs within cutoff will be returned.\n"
 396                 << "  Actual: Position " << j->refIndex << " is not found.";
 397                 break;
 398             }
 399         }
 400     }
 401 }
 402
 403 /********************************************************************
 404  * Test data generation
 405  */
 406
 407 class TrivialTestData
 408 {
 409     public:
 410         static const NeighborhoodSearchTestData &get()
 411         {
 412             static TrivialTestData singleton;
 413             return singleton.data_;
 414         }
 415
 416         TrivialTestData() : data_(12345, 1.0)
 417         {
 418             data_.box_[XX][XX] = 5.0;
 419             data_.box_[YY][YY] = 5.0;
 420             data_.box_[ZZ][ZZ] = 5.0;
 421             data_.generateRandomRefPositions(10);
 422             data_.generateRandomTestPositions(5);
 423             set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
 424             data_.computeReferences(&data_.pbc_);
 425         }
 426
 427     private:
 428         NeighborhoodSearchTestData data_;
 429 };
 430
 431 class RandomBoxFullPBCData
 432 {
 433     public:
 434         static const NeighborhoodSearchTestData &get()
 435         {
 436             static RandomBoxFullPBCData singleton;
 437             return singleton.data_;
 438         }
 439
 440         RandomBoxFullPBCData() : data_(12345, 1.0)
 441         {
 442             data_.box_[XX][XX] = 10.0;
 443             data_.box_[YY][YY] = 5.0;
 444             data_.box_[ZZ][ZZ] = 7.0;
 445             // TODO: Consider whether manually picking some positions would give better
 446             // test coverage.
 447             data_.generateRandomRefPositions(1000);
 448             data_.generateRandomTestPositions(100);
 449             set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
 450             data_.computeReferences(&data_.pbc_);
 451         }
 452
 453     private:
 454         NeighborhoodSearchTestData data_;
 455 };
 456
 457 class RandomBoxXYFullPBCData
 458 {
 459     public:
 460         static const NeighborhoodSearchTestData &get()
 461         {
 462             static RandomBoxXYFullPBCData singleton;
 463             return singleton.data_;
 464         }
 465
 466         RandomBoxXYFullPBCData() : data_(54321, 1.0)
 467         {
 468             data_.box_[XX][XX] = 10.0;
 469             data_.box_[YY][YY] = 5.0;
 470             data_.box_[ZZ][ZZ] = 7.0;
 471             // TODO: Consider whether manually picking some positions would give better
 472             // test coverage.
 473             data_.generateRandomRefPositions(1000);
 474             data_.generateRandomTestPositions(100);
 475             set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
 476             data_.computeReferencesXY(&data_.pbc_);
 477         }
 478
 479     private:
 480         NeighborhoodSearchTestData data_;
 481 };
 482
 483 class RandomTriclinicFullPBCData
 484 {
 485     public:
 486         static const NeighborhoodSearchTestData &get()
 487         {
 488             static RandomTriclinicFullPBCData singleton;
 489             return singleton.data_;
 490         }
 491
 492         RandomTriclinicFullPBCData() : data_(12345, 1.0)
 493         {
 494             data_.box_[XX][XX] = 5.0;
 495             data_.box_[YY][XX] = 2.5;
 496             data_.box_[YY][YY] = 2.5*sqrt(3.0);
 497             data_.box_[ZZ][XX] = 2.5;
 498             data_.box_[ZZ][YY] = 2.5*sqrt(1.0/3.0);
 499             data_.box_[ZZ][ZZ] = 5.0*sqrt(2.0/3.0);
 500             // TODO: Consider whether manually picking some positions would give better
 501             // test coverage.
 502             data_.generateRandomRefPositions(1000);
 503             data_.generateRandomTestPositions(100);
 504             set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
 505             data_.computeReferences(&data_.pbc_);
 506         }
 507
 508     private:
 509         NeighborhoodSearchTestData data_;
 510 };
 511
 512 class RandomBox2DPBCData
 513 {
 514     public:
 515         static const NeighborhoodSearchTestData &get()
 516         {
 517             static RandomBox2DPBCData singleton;
 518             return singleton.data_;
 519         }
 520
 521         RandomBox2DPBCData() : data_(12345, 1.0)
 522         {
 523             data_.box_[XX][XX] = 10.0;
 524             data_.box_[YY][YY] = 7.0;
 525             data_.box_[ZZ][ZZ] = 5.0;
 526             // TODO: Consider whether manually picking some positions would give better
 527             // test coverage.
 528             data_.generateRandomRefPositions(1000);
 529             data_.generateRandomTestPositions(100);
 530             set_pbc(&data_.pbc_, epbcXY, data_.box_);
 531             data_.computeReferences(&data_.pbc_);
 532         }
 533
 534     private:
 535         NeighborhoodSearchTestData data_;
 536 };
 537
 538 /********************************************************************
 539  * Actual tests
 540  */
 541
 542 TEST_F(NeighborhoodSearchTest, SimpleSearch)
 543 {
 544     const NeighborhoodSearchTestData &data = RandomBoxFullPBCData::get();
 545
 546     nb_.setCutoff(data.cutoff_);
 547     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Simple);
 548     gmx::AnalysisNeighborhoodSearch search =
 549         nb_.initSearch(&data.pbc_, data.refPositions());
 550     ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Simple, search.mode());
 551
 552     testIsWithin(&search, data);
 553     testMinimumDistance(&search, data);
 554     testNearestPoint(&search, data);
 555     testPairSearch(&search, data);
 556 }
 557
 558 TEST_F(NeighborhoodSearchTest, GridSearchBox)
 559 {
 560     const NeighborhoodSearchTestData &data = RandomBoxFullPBCData::get();
 561
 562     nb_.setCutoff(data.cutoff_);
 563     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
 564     gmx::AnalysisNeighborhoodSearch search =
 565         nb_.initSearch(&data.pbc_, data.refPositions());
 566     ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
 567
 568     testIsWithin(&search, data);
 569     testMinimumDistance(&search, data);
 570     testNearestPoint(&search, data);
 571     testPairSearch(&search, data);
 572 }
 573
 574 TEST_F(NeighborhoodSearchTest, GridSearchTriclinic)
 575 {
 576     const NeighborhoodSearchTestData &data = RandomTriclinicFullPBCData::get();
 577
 578     nb_.setCutoff(data.cutoff_);
 579     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
 580     gmx::AnalysisNeighborhoodSearch search =
 581         nb_.initSearch(&data.pbc_, data.refPositions());
 582     ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
 583
 584     testPairSearch(&search, data);
 585 }
 586
 587 TEST_F(NeighborhoodSearchTest, GridSearch2DPBC)
 588 {
 589     const NeighborhoodSearchTestData &data = RandomBox2DPBCData::get();
 590
 591     nb_.setCutoff(data.cutoff_);
 592     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
 593     gmx::AnalysisNeighborhoodSearch search =
 594         nb_.initSearch(&data.pbc_, data.refPositions());
 595     // Currently, grid searching not supported with 2D PBC.
 596     //ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
 597
 598     testIsWithin(&search, data);
 599     testMinimumDistance(&search, data);
 600     testNearestPoint(&search, data);
 601     testPairSearch(&search, data);
 602 }
 603
 604 TEST_F(NeighborhoodSearchTest, GridSearchXYBox)
 605 {
 606     const NeighborhoodSearchTestData &data = RandomBoxXYFullPBCData::get();
 607
 608     nb_.setCutoff(data.cutoff_);
 609     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
 610     nb_.setXYMode(true);
 611     gmx::AnalysisNeighborhoodSearch search =
 612         nb_.initSearch(&data.pbc_, data.refPositions());
 613     // Currently, grid searching not supported with XY.
 614     //ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
 615
 616     testIsWithin(&search, data);
 617     testMinimumDistance(&search, data);
 618     testNearestPoint(&search, data);
 619     testPairSearch(&search, data);
 620 }
 621
 622 TEST_F(NeighborhoodSearchTest, HandlesConcurrentSearches)
 623 {
 624     const NeighborhoodSearchTestData &data = TrivialTestData::get();
 625
 626     nb_.setCutoff(data.cutoff_);
 627     gmx::AnalysisNeighborhoodSearch search1 =
 628         nb_.initSearch(&data.pbc_, data.refPositions());
 629     gmx::AnalysisNeighborhoodSearch search2 =
 630         nb_.initSearch(&data.pbc_, data.refPositions());
 631
 632     gmx::AnalysisNeighborhoodPairSearch pairSearch1 =
 633         search1.startPairSearch(data.testPosition(0));
 634     gmx::AnalysisNeighborhoodPairSearch pairSearch2 =
 635         search1.startPairSearch(data.testPosition(1));
 636
 637     testPairSearch(&search2, data);
 638
 639     gmx::AnalysisNeighborhoodPair pair;
 640     pairSearch1.findNextPair(&pair);
 641     EXPECT_EQ(0, pair.testIndex());
 642     {
 643         NeighborhoodSearchTestData::RefPair searchPair(pair.refIndex(), sqrt(pair.distance2()));
 644         EXPECT_TRUE(data.containsPair(0, searchPair));
 645     }
 646
 647     pairSearch2.findNextPair(&pair);
 648     EXPECT_EQ(1, pair.testIndex());
 649     {
 650         NeighborhoodSearchTestData::RefPair searchPair(pair.refIndex(), sqrt(pair.distance2()));
 651         EXPECT_TRUE(data.containsPair(1, searchPair));
 652     }
 653 }
 654
 655 TEST_F(NeighborhoodSearchTest, HandlesSkippingPairs)
 656 {
 657     const NeighborhoodSearchTestData &data = TrivialTestData::get();
 658
 659     nb_.setCutoff(data.cutoff_);
 660     gmx::AnalysisNeighborhoodSearch     search =
 661         nb_.initSearch(&data.pbc_, data.refPositions());
 662     gmx::AnalysisNeighborhoodPairSearch pairSearch =
 663         search.startPairSearch(data.testPositions());
 664     gmx::AnalysisNeighborhoodPair       pair;
 665     // TODO: This test needs to be adjusted if the grid search gets optimized
 666     // to loop over the test positions in cell order (first, the ordering
 667     // assumption here breaks, and second, it then needs to be tested
 668     // separately for simple and grid searches).
 669     int currentIndex = 0;
 670     while (pairSearch.findNextPair(&pair))
 671     {
 672         while (currentIndex < pair.testIndex())
 673         {
 674             ++currentIndex;
 675         }
 676         EXPECT_EQ(currentIndex, pair.testIndex());
 677         NeighborhoodSearchTestData::RefPair searchPair(pair.refIndex(), sqrt(pair.distance2()));
 678         EXPECT_TRUE(data.containsPair(currentIndex, searchPair));
 679         pairSearch.skipRemainingPairsForTestPosition();
 680         ++currentIndex;
 681     }
 682 }
 683
 684 } // namespace