2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2013, by the GROMACS development team, led by
5 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
6 * others, as listed in the AUTHORS file in the top-level source
7 * directory and at http://www.gromacs.org.
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.
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.
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.
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.
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.
37 * Tests selection neighborhood searching.
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.
44 * \author Teemu Murtola <teemu.murtola@gmail.com>
45 * \ingroup module_selection
47 #include <gtest/gtest.h>
55 #include "gromacs/legacyheaders/gmx_random.h"
56 #include "gromacs/legacyheaders/pbc.h"
57 #include "gromacs/legacyheaders/smalloc.h"
58 #include "gromacs/legacyheaders/vec.h"
60 #include "gromacs/selection/nbsearch.h"
62 #include "testutils/testasserts.h"
67 /********************************************************************
68 * NeighborhoodSearchTestData
71 class NeighborhoodSearchTestData
76 TestPosition() : refMinDist(0.0), refNearestPoint(-1)
80 explicit TestPosition(const rvec x)
81 : refMinDist(0.0), refNearestPoint(-1)
83 copy_rvec(x, this->x);
89 std::set<int> refPairs;
91 typedef std::vector<TestPosition> TestPositionList;
93 NeighborhoodSearchTestData(int seed, real cutoff);
94 ~NeighborhoodSearchTestData();
96 void addTestPosition(const rvec x)
98 testPositions_.push_back(TestPosition(x));
100 void generateRandomPosition(rvec x);
101 void generateRandomRefPositions(int count);
102 void generateRandomTestPositions(int count);
103 void computeReferences(t_pbc *pbc);
111 TestPositionList testPositions_;
114 NeighborhoodSearchTestData::NeighborhoodSearchTestData(int seed, real cutoff)
115 : rng_(NULL), cutoff_(cutoff), refPosCount_(0), refPos_(NULL)
117 // TODO: Handle errors.
118 rng_ = gmx_rng_init(seed);
120 set_pbc(&pbc_, epbcNONE, box_);
123 NeighborhoodSearchTestData::~NeighborhoodSearchTestData()
127 gmx_rng_destroy(rng_);
132 void NeighborhoodSearchTestData::generateRandomPosition(rvec x)
135 fx[XX] = gmx_rng_uniform_real(rng_);
136 fx[YY] = gmx_rng_uniform_real(rng_);
137 fx[ZZ] = gmx_rng_uniform_real(rng_);
139 // Add a small displacement to allow positions outside the box
140 x[XX] += 0.2 * gmx_rng_uniform_real(rng_) - 0.1;
141 x[YY] += 0.2 * gmx_rng_uniform_real(rng_) - 0.1;
142 x[ZZ] += 0.2 * gmx_rng_uniform_real(rng_) - 0.1;
145 void NeighborhoodSearchTestData::generateRandomRefPositions(int count)
147 refPosCount_ = count;
148 snew(refPos_, refPosCount_);
149 for (int i = 0; i < refPosCount_; ++i)
151 generateRandomPosition(refPos_[i]);
155 void NeighborhoodSearchTestData::generateRandomTestPositions(int count)
157 testPositions_.reserve(count);
158 for (int i = 0; i < count; ++i)
161 generateRandomPosition(x);
166 void NeighborhoodSearchTestData::computeReferences(t_pbc *pbc)
168 real cutoff = cutoff_;
171 cutoff = std::numeric_limits<real>::max();
173 TestPositionList::iterator i;
174 for (i = testPositions_.begin(); i != testPositions_.end(); ++i)
176 i->refMinDist = cutoff;
177 i->refNearestPoint = -1;
179 for (int j = 0; j < refPosCount_; ++j)
184 pbc_dx(pbc, i->x, refPos_[j], dx);
188 rvec_sub(i->x, refPos_[j], dx);
190 const real dist = norm(dx);
191 if (dist < i->refMinDist)
193 i->refMinDist = dist;
194 i->refNearestPoint = j;
198 i->refPairs.insert(j);
204 /********************************************************************
205 * NeighborhoodSearchTest
208 class NeighborhoodSearchTest : public ::testing::Test
211 void testIsWithin(gmx::AnalysisNeighborhoodSearch *search,
212 const NeighborhoodSearchTestData &data);
213 void testMinimumDistance(gmx::AnalysisNeighborhoodSearch *search,
214 const NeighborhoodSearchTestData &data);
215 void testNearestPoint(gmx::AnalysisNeighborhoodSearch *search,
216 const NeighborhoodSearchTestData &data);
217 void testPairSearch(gmx::AnalysisNeighborhoodSearch *search,
218 const NeighborhoodSearchTestData &data);
220 gmx::AnalysisNeighborhood nb_;
223 void NeighborhoodSearchTest::testIsWithin(
224 gmx::AnalysisNeighborhoodSearch *search,
225 const NeighborhoodSearchTestData &data)
227 NeighborhoodSearchTestData::TestPositionList::const_iterator i;
228 for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
230 const bool bWithin = (i->refMinDist <= data.cutoff_);
231 EXPECT_EQ(bWithin, search->isWithin(i->x))
232 << "Distance is " << i->refMinDist;
236 void NeighborhoodSearchTest::testMinimumDistance(
237 gmx::AnalysisNeighborhoodSearch *search,
238 const NeighborhoodSearchTestData &data)
240 NeighborhoodSearchTestData::TestPositionList::const_iterator i;
241 for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
243 const real refDist = i->refMinDist;
244 EXPECT_NEAR_REL(refDist, search->minimumDistance(i->x), 20*GMX_REAL_EPS);
248 void NeighborhoodSearchTest::testNearestPoint(
249 gmx::AnalysisNeighborhoodSearch *search,
250 const NeighborhoodSearchTestData &data)
252 NeighborhoodSearchTestData::TestPositionList::const_iterator i;
253 for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
255 const gmx::AnalysisNeighborhoodPair pair = search->nearestPoint(i->x);
258 EXPECT_EQ(i->refNearestPoint, pair.refIndex());
259 EXPECT_EQ(0, pair.testIndex());
263 EXPECT_EQ(i->refNearestPoint, -1);
268 void NeighborhoodSearchTest::testPairSearch(
269 gmx::AnalysisNeighborhoodSearch *search,
270 const NeighborhoodSearchTestData &data)
272 NeighborhoodSearchTestData::TestPositionList::const_iterator i;
273 for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
275 std::set<int> checkSet = i->refPairs;
276 gmx::AnalysisNeighborhoodPairSearch pairSearch =
277 search->startPairSearch(i->x);
278 gmx::AnalysisNeighborhoodPair pair;
279 while (pairSearch.findNextPair(&pair))
281 EXPECT_EQ(0, pair.testIndex());
282 if (checkSet.erase(pair.refIndex()) == 0)
284 // TODO: Check whether the same pair was returned more than
285 // once and give a better error message if so.
287 << "Expected: Position " << pair.refIndex()
288 << " is within cutoff.\n"
289 << " Actual: It is not.";
292 EXPECT_TRUE(checkSet.empty()) << "Some positions were not returned by the pair search.";
296 /********************************************************************
297 * Test data generation
300 class TrivialTestData
303 static const NeighborhoodSearchTestData &get()
305 static TrivialTestData singleton;
306 return singleton.data_;
309 TrivialTestData() : data_(12345, 1.0)
311 data_.box_[XX][XX] = 5.0;
312 data_.box_[YY][YY] = 5.0;
313 data_.box_[ZZ][ZZ] = 5.0;
314 data_.generateRandomRefPositions(10);
315 data_.generateRandomTestPositions(5);
316 set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
317 data_.computeReferences(&data_.pbc_);
321 NeighborhoodSearchTestData data_;
324 class RandomBoxFullPBCData
327 static const NeighborhoodSearchTestData &get()
329 static RandomBoxFullPBCData singleton;
330 return singleton.data_;
333 RandomBoxFullPBCData() : data_(12345, 1.0)
335 data_.box_[XX][XX] = 10.0;
336 data_.box_[YY][YY] = 5.0;
337 data_.box_[ZZ][ZZ] = 7.0;
338 // TODO: Consider whether manually picking some positions would give better
340 data_.generateRandomRefPositions(1000);
341 data_.generateRandomTestPositions(100);
342 set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
343 data_.computeReferences(&data_.pbc_);
347 NeighborhoodSearchTestData data_;
350 class RandomTriclinicFullPBCData
353 static const NeighborhoodSearchTestData &get()
355 static RandomTriclinicFullPBCData singleton;
356 return singleton.data_;
359 RandomTriclinicFullPBCData() : data_(12345, 1.0)
361 data_.box_[XX][XX] = 5.0;
362 data_.box_[YY][XX] = 2.5;
363 data_.box_[YY][YY] = 2.5*sqrt(3.0);
364 data_.box_[ZZ][XX] = 2.5;
365 data_.box_[ZZ][YY] = 2.5*sqrt(1.0/3.0);
366 data_.box_[ZZ][ZZ] = 5.0*sqrt(2.0/3.0);
367 // TODO: Consider whether manually picking some positions would give better
369 data_.generateRandomRefPositions(1000);
370 data_.generateRandomTestPositions(100);
371 set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
372 data_.computeReferences(&data_.pbc_);
376 NeighborhoodSearchTestData data_;
379 class RandomBox2DPBCData
382 static const NeighborhoodSearchTestData &get()
384 static RandomBox2DPBCData singleton;
385 return singleton.data_;
388 RandomBox2DPBCData() : data_(12345, 1.0)
390 data_.box_[XX][XX] = 10.0;
391 data_.box_[YY][YY] = 7.0;
392 data_.box_[ZZ][ZZ] = 5.0;
393 // TODO: Consider whether manually picking some positions would give better
395 data_.generateRandomRefPositions(1000);
396 data_.generateRandomTestPositions(100);
397 set_pbc(&data_.pbc_, epbcXY, data_.box_);
398 data_.computeReferences(&data_.pbc_);
402 NeighborhoodSearchTestData data_;
405 /********************************************************************
409 TEST_F(NeighborhoodSearchTest, SimpleSearch)
411 const NeighborhoodSearchTestData &data = RandomBoxFullPBCData::get();
413 nb_.setCutoff(data.cutoff_);
414 nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Simple);
415 gmx::AnalysisNeighborhoodSearch search =
416 nb_.initSearch(&data.pbc_, data.refPosCount_, data.refPos_);
417 ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Simple, search.mode());
419 testIsWithin(&search, data);
420 testMinimumDistance(&search, data);
421 testNearestPoint(&search, data);
422 testPairSearch(&search, data);
425 TEST_F(NeighborhoodSearchTest, GridSearchBox)
427 const NeighborhoodSearchTestData &data = RandomBoxFullPBCData::get();
429 nb_.setCutoff(data.cutoff_);
430 nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
431 gmx::AnalysisNeighborhoodSearch search =
432 nb_.initSearch(&data.pbc_, data.refPosCount_, data.refPos_);
433 ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
435 testIsWithin(&search, data);
436 testMinimumDistance(&search, data);
437 testNearestPoint(&search, data);
438 testPairSearch(&search, data);
441 TEST_F(NeighborhoodSearchTest, GridSearchTriclinic)
443 const NeighborhoodSearchTestData &data = RandomTriclinicFullPBCData::get();
445 nb_.setCutoff(data.cutoff_);
446 nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
447 gmx::AnalysisNeighborhoodSearch search =
448 nb_.initSearch(&data.pbc_, data.refPosCount_, data.refPos_);
449 ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
451 testPairSearch(&search, data);
454 TEST_F(NeighborhoodSearchTest, GridSearch2DPBC)
456 const NeighborhoodSearchTestData &data = RandomBox2DPBCData::get();
458 nb_.setCutoff(data.cutoff_);
459 nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
460 gmx::AnalysisNeighborhoodSearch search =
461 nb_.initSearch(&data.pbc_, data.refPosCount_, data.refPos_);
462 // Currently, grid searching not supported with 2D PBC.
463 //ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
465 testIsWithin(&search, data);
466 testMinimumDistance(&search, data);
467 testNearestPoint(&search, data);
468 testPairSearch(&search, data);
471 TEST_F(NeighborhoodSearchTest, HandlesConcurrentSearches)
473 const NeighborhoodSearchTestData &data = TrivialTestData::get();
475 nb_.setCutoff(data.cutoff_);
476 gmx::AnalysisNeighborhoodSearch search1 =
477 nb_.initSearch(&data.pbc_, data.refPosCount_, data.refPos_);
478 gmx::AnalysisNeighborhoodSearch search2 =
479 nb_.initSearch(&data.pbc_, data.refPosCount_, data.refPos_);
481 gmx::AnalysisNeighborhoodPairSearch pairSearch1 =
482 search1.startPairSearch(data.testPositions_[0].x);
483 gmx::AnalysisNeighborhoodPairSearch pairSearch2 =
484 search1.startPairSearch(data.testPositions_[1].x);
486 testPairSearch(&search2, data);
488 gmx::AnalysisNeighborhoodPair pair;
489 pairSearch1.findNextPair(&pair);
490 EXPECT_EQ(0, pair.testIndex());
491 EXPECT_TRUE(data.testPositions_[0].refPairs.count(pair.refIndex()) == 1);
493 pairSearch2.findNextPair(&pair);
494 EXPECT_EQ(0, pair.testIndex());
495 EXPECT_TRUE(data.testPositions_[1].refPairs.count(pair.refIndex()) == 1);