src/gromacs/selection/tests/poscalc.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2013, 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 the position mapping engine.
  38  *
  39  * \author Teemu Murtola <teemu.murtola@gmail.com>
  40  * \ingroup module_selection
  41  */
  42 #include <gtest/gtest.h>
  43
  44 #include <vector>
  45
  46 #include "gromacs/legacyheaders/smalloc.h"
  47 #include "gromacs/legacyheaders/typedefs.h"
  48 #include "gromacs/legacyheaders/vec.h"
  49
  50 #include "gromacs/selection/indexutil.h"
  51 #include "gromacs/selection/poscalc.h"
  52 #include "gromacs/selection/position.h"
  53 #include "gromacs/utility/uniqueptr.h"
  54
  55 #include "testutils/refdata.h"
  56
  57 #include "toputils.h"
  58
  59 namespace
  60 {
  61
  62 /********************************************************************
  63  * PositionCalculationTest
  64  */
  65
  66 class PositionCalculationTest : public ::testing::Test
  67 {
  68     public:
  69         PositionCalculationTest();
  70         ~PositionCalculationTest();
  71
  72         void generateCoordinates();
  73
  74         gmx_ana_poscalc_t *createCalculation(e_poscalc_t type, int flags);
  75         void setMaximumGroup(gmx_ana_poscalc_t *pc,
  76                              int count, const int atoms[]);
  77         gmx_ana_pos_t *initPositions(gmx_ana_poscalc_t *pc, const char *name);
  78
  79         void checkInitialized();
  80         void updateAndCheck(gmx_ana_poscalc_t *pc, gmx_ana_pos_t *p,
  81                             int count, const int atoms[],
  82                             gmx::test::TestReferenceChecker *checker,
  83                             const char *name);
  84
  85         void testSingleStatic(e_poscalc_t type, int flags, bool bExpectTop,
  86                               int atomCount, const int atoms[]);
  87         void testSingleDynamic(e_poscalc_t type, int flags, bool bExpectTop,
  88                                int initCount, const int initAtoms[],
  89                                int evalCount, const int evalAtoms[]);
  90
  91         template <int count>
  92         void setMaximumGroup(gmx_ana_poscalc_t *pc, const int (&atoms)[count])
  93         {
  94             setMaximumGroup(pc, count, atoms);
  95         }
  96         template <int count>
  97         void updateAndCheck(gmx_ana_poscalc_t *pc, gmx_ana_pos_t *p,
  98                             const int (&atoms)[count],
  99                             gmx::test::TestReferenceChecker *checker,
 100                             const char *name)
 101         {
 102             updateAndCheck(pc, p, count, atoms, checker, name);
 103         }
 104         template <int atomCount>
 105         void testSingleStatic(e_poscalc_t type, int flags, bool bExpectTop,
 106                               const int (&atoms)[atomCount])
 107         {
 108             testSingleStatic(type, flags, bExpectTop, atomCount, atoms);
 109         }
 110         template <int initCount, int evalCount>
 111         void testSingleDynamic(e_poscalc_t type, int flags, bool bExpectTop,
 112                                const int (&initAtoms)[initCount],
 113                                const int (&evalAtoms)[evalCount])
 114         {
 115             testSingleDynamic(type, flags, bExpectTop,
 116                               initCount, initAtoms, evalCount, evalAtoms);
 117         }
 118
 119         gmx::test::TestReferenceData        data_;
 120         gmx::test::TestReferenceChecker     checker_;
 121         gmx::test::TopologyManager          topManager_;
 122         gmx::PositionCalculationCollection  pcc_;
 123
 124     private:
 125         typedef gmx::gmx_unique_ptr<gmx_ana_pos_t>::type PositionPointer;
 126
 127         struct PositionTest
 128         {
 129             PositionTest(PositionPointer pos, gmx_ana_poscalc_t *pc,
 130                          const char *name)
 131                 : pos(gmx::move(pos)), pc(pc), name(name)
 132             {
 133             }
 134
 135             PositionPointer                 pos;
 136             gmx_ana_poscalc_t              *pc;
 137             const char                     *name;
 138         };
 139
 140         typedef std::vector<PositionTest> PositionTestList;
 141
 142         void setTopologyIfRequired();
 143         void checkPositions(gmx::test::TestReferenceChecker *checker,
 144                             const char *name, gmx_ana_pos_t *p,
 145                             bool bCoordinates);
 146
 147         std::vector<gmx_ana_poscalc_t *>    pcList_;
 148         PositionTestList                    posList_;
 149         bool                                bTopSet_;
 150 };
 151
 152 PositionCalculationTest::PositionCalculationTest()
 153     : checker_(data_.rootChecker()), bTopSet_(false)
 154 {
 155     topManager_.requestFrame();
 156 }
 157
 158 PositionCalculationTest::~PositionCalculationTest()
 159 {
 160     std::vector<gmx_ana_poscalc_t *>::reverse_iterator pci;
 161     for (pci = pcList_.rbegin(); pci != pcList_.rend(); ++pci)
 162     {
 163         gmx_ana_poscalc_free(*pci);
 164     }
 165 }
 166
 167 void PositionCalculationTest::generateCoordinates()
 168 {
 169     t_topology *top   = topManager_.topology();
 170     t_trxframe *frame = topManager_.frame();
 171     for (int i = 0; i < top->atoms.nr; ++i)
 172     {
 173         frame->x[i][XX] = i;
 174         frame->x[i][YY] = top->atoms.atom[i].resind;
 175         frame->x[i][ZZ] = 0.0;
 176         if (frame->bV)
 177         {
 178             copy_rvec(frame->x[i], frame->v[i]);
 179             frame->v[i][ZZ] = 1.0;
 180         }
 181         if (frame->bF)
 182         {
 183             copy_rvec(frame->x[i], frame->f[i]);
 184             frame->f[i][ZZ] = -1.0;
 185         }
 186     }
 187 }
 188
 189 gmx_ana_poscalc_t *
 190 PositionCalculationTest::createCalculation(e_poscalc_t type, int flags)
 191 {
 192     pcList_.reserve(pcList_.size() + 1);
 193     pcList_.push_back(pcc_.createCalculation(type, flags));
 194     return pcList_.back();
 195 }
 196
 197 void PositionCalculationTest::setMaximumGroup(gmx_ana_poscalc_t *pc,
 198                                               int count, const int atoms[])
 199 {
 200     setTopologyIfRequired();
 201     gmx_ana_index_t g;
 202     g.isize = count;
 203     g.index = const_cast<int *>(atoms);
 204     gmx_ana_poscalc_set_maxindex(pc, &g);
 205 }
 206
 207 gmx_ana_pos_t *
 208 PositionCalculationTest::initPositions(gmx_ana_poscalc_t *pc, const char *name)
 209 {
 210     posList_.reserve(posList_.size() + 1);
 211     PositionPointer p(new gmx_ana_pos_t());
 212     gmx_ana_pos_t  *result = p.get();
 213     posList_.push_back(PositionTest(gmx::move(p), pc, name));
 214     gmx_ana_poscalc_init_pos(pc, result);
 215     return result;
 216 }
 217
 218 void PositionCalculationTest::checkInitialized()
 219 {
 220     gmx::test::TestReferenceChecker  compound(
 221             checker_.checkCompound("InitializedPositions", NULL));
 222     PositionTestList::const_iterator pi;
 223     for (pi = posList_.begin(); pi != posList_.end(); ++pi)
 224     {
 225         checkPositions(&compound, pi->name, pi->pos.get(), false);
 226     }
 227 }
 228
 229 void PositionCalculationTest::updateAndCheck(
 230         gmx_ana_poscalc_t *pc, gmx_ana_pos_t *p, int count, const int atoms[],
 231         gmx::test::TestReferenceChecker *checker, const char *name)
 232 {
 233     gmx_ana_index_t g;
 234     g.isize = count;
 235     g.index = const_cast<int *>(atoms);
 236     gmx_ana_poscalc_update(pc, p, &g, topManager_.frame(), NULL);
 237     checkPositions(checker, name, p, true);
 238 }
 239
 240 void PositionCalculationTest::testSingleStatic(
 241         e_poscalc_t type, int flags, bool bExpectTop,
 242         int atomCount, const int atoms[])
 243 {
 244     t_trxframe *frame = topManager_.frame();
 245     if (frame->bV)
 246     {
 247         flags |= POS_VELOCITIES;
 248     }
 249     if (frame->bF)
 250     {
 251         flags |= POS_FORCES;
 252     }
 253     gmx_ana_poscalc_t *pc = createCalculation(type, flags);
 254     EXPECT_EQ(bExpectTop, gmx_ana_poscalc_requires_top(pc));
 255     setMaximumGroup(pc, atomCount, atoms);
 256     gmx_ana_pos_t *p = initPositions(pc, NULL);
 257     checkInitialized();
 258     {
 259         pcc_.initEvaluation();
 260         pcc_.initFrame();
 261         generateCoordinates();
 262         gmx::test::TestReferenceChecker frameCompound(
 263                 checker_.checkCompound("EvaluatedPositions", "Frame0"));
 264         updateAndCheck(pc, p, atomCount, atoms, &frameCompound, NULL);
 265     }
 266 }
 267
 268 void PositionCalculationTest::testSingleDynamic(
 269         e_poscalc_t type, int flags, bool bExpectTop,
 270         int initCount, const int initAtoms[],
 271         int evalCount, const int evalAtoms[])
 272 {
 273     gmx_ana_poscalc_t *pc = createCalculation(type, flags | POS_DYNAMIC);
 274     EXPECT_EQ(bExpectTop, gmx_ana_poscalc_requires_top(pc));
 275     setMaximumGroup(pc, initCount, initAtoms);
 276     gmx_ana_pos_t *p = initPositions(pc, NULL);
 277     checkInitialized();
 278     {
 279         pcc_.initEvaluation();
 280         pcc_.initFrame();
 281         generateCoordinates();
 282         gmx::test::TestReferenceChecker frameCompound(
 283                 checker_.checkCompound("EvaluatedPositions", "Frame0"));
 284         updateAndCheck(pc, p, evalCount, evalAtoms, &frameCompound, NULL);
 285     }
 286 }
 287
 288 void PositionCalculationTest::setTopologyIfRequired()
 289 {
 290     if (bTopSet_)
 291     {
 292         return;
 293     }
 294     std::vector<gmx_ana_poscalc_t *>::const_iterator pci;
 295     for (pci = pcList_.begin(); pci != pcList_.end(); ++pci)
 296     {
 297         if (gmx_ana_poscalc_requires_top(*pci))
 298         {
 299             bTopSet_ = true;
 300             pcc_.setTopology(topManager_.topology());
 301             return;
 302         }
 303     }
 304 }
 305
 306 void PositionCalculationTest::checkPositions(
 307         gmx::test::TestReferenceChecker *checker,
 308         const char *name, gmx_ana_pos_t *p, bool bCoordinates)
 309 {
 310     gmx::test::TestReferenceChecker compound(
 311             checker->checkCompound("Positions", name));
 312     compound.checkInteger(p->count(), "Count");
 313     const char *type = "???";
 314     switch (p->m.type)
 315     {
 316         case INDEX_UNKNOWN: type = "unknown";   break;
 317         case INDEX_ATOM:    type = "atoms";     break;
 318         case INDEX_RES:     type = "residues";  break;
 319         case INDEX_MOL:     type = "molecules"; break;
 320         case INDEX_ALL:     type = "single";    break;
 321     }
 322     compound.checkString(type, "Type");
 323     compound.checkSequenceArray(p->count() + 1, p->m.mapb.index, "Block");
 324     for (int i = 0; i < p->count(); ++i)
 325     {
 326         gmx::test::TestReferenceChecker posCompound(
 327                 compound.checkCompound("Position", NULL));
 328         posCompound.checkSequence(&p->m.mapb.a[p->m.mapb.index[i]],
 329                                   &p->m.mapb.a[p->m.mapb.index[i+1]],
 330                                   "Atoms");
 331         posCompound.checkInteger(p->m.refid[i], "RefId");
 332         if (bCoordinates)
 333         {
 334             posCompound.checkVector(p->x[i], "Coordinates");
 335         }
 336         if (bCoordinates && p->v != NULL)
 337         {
 338             posCompound.checkVector(p->v[i], "Velocity");
 339         }
 340         if (bCoordinates && p->f != NULL)
 341         {
 342             posCompound.checkVector(p->f[i], "Force");
 343         }
 344         int originalIdIndex = (p->m.refid[i] != -1 ? p->m.refid[i] : i);
 345         EXPECT_EQ(p->m.orgid[originalIdIndex], p->m.mapid[i]);
 346     }
 347 }
 348
 349 /********************************************************************
 350  * Actual tests
 351  */
 352
 353 TEST_F(PositionCalculationTest, ComputesAtomPositions)
 354 {
 355     const int group[] = { 0, 1, 2, 3 };
 356     topManager_.requestVelocities();
 357     topManager_.requestForces();
 358     topManager_.initAtoms(4);
 359     testSingleStatic(POS_ATOM, 0, false, group);
 360 }
 361
 362 TEST_F(PositionCalculationTest, ComputesResidueCOGPositions)
 363 {
 364     const int group[] = { 0, 1, 2, 3, 4, 8 };
 365     topManager_.requestVelocities();
 366     topManager_.requestForces();
 367     topManager_.initAtoms(9);
 368     topManager_.initUniformResidues(3);
 369     testSingleStatic(POS_RES, 0, true, group);
 370 }
 371
 372 TEST_F(PositionCalculationTest, ComputesResidueCOMPositions)
 373 {
 374     const int group[] = { 0, 1, 2, 3, 4, 8 };
 375     topManager_.requestVelocities();
 376     topManager_.requestForces();
 377     topManager_.initAtoms(9);
 378     topManager_.initUniformResidues(3);
 379     testSingleStatic(POS_RES, POS_MASS, true, group);
 380 }
 381
 382 TEST_F(PositionCalculationTest, ComputesGroupCOGPositions)
 383 {
 384     const int group[] = { 0, 1, 2, 3, 4, 5, 6, 7, 8 };
 385     topManager_.requestVelocities();
 386     topManager_.requestForces();
 387     topManager_.initAtoms(9);
 388     // Topology (masses) is requires for computing the force
 389     testSingleStatic(POS_ALL, 0, true, group);
 390 }
 391
 392 TEST_F(PositionCalculationTest, ComputesGroupCOMPositions)
 393 {
 394     const int group[] = { 0, 1, 2, 3, 4, 5, 6, 7, 8 };
 395     topManager_.requestVelocities();
 396     topManager_.requestForces();
 397     topManager_.initAtoms(9);
 398     testSingleStatic(POS_ALL, POS_MASS, true, group);
 399 }
 400
 401 TEST_F(PositionCalculationTest, ComputesPositionsWithCompleteWhole)
 402 {
 403     const int group[] = { 0, 1, 2, 3, 4, 8 };
 404     topManager_.initAtoms(9);
 405     topManager_.initUniformResidues(3);
 406     testSingleStatic(POS_RES, POS_COMPLWHOLE, true, group);
 407 }
 408
 409 TEST_F(PositionCalculationTest, ComputesPositionsWithCompleteMax)
 410 {
 411     const int maxGroup[]  = { 0, 1, 4, 5, 6, 8 };
 412     const int evalGroup[] = { 0, 1, 5, 6 };
 413     topManager_.initAtoms(9);
 414     topManager_.initUniformResidues(3);
 415     testSingleDynamic(POS_RES, POS_COMPLMAX, true, maxGroup, evalGroup);
 416 }
 417
 418 TEST_F(PositionCalculationTest, ComputesPositionMask)
 419 {
 420     const int maxGroup[]  = { 0, 1, 2, 3, 4, 5 };
 421     const int evalGroup[] = { 1, 2, 4 };
 422     topManager_.initAtoms(6);
 423     testSingleDynamic(POS_ATOM, POS_MASKONLY, false, maxGroup, evalGroup);
 424 }
 425
 426 // TODO: Check for POS_ALL_PBC
 427
 428 TEST_F(PositionCalculationTest, HandlesIdenticalStaticCalculations)
 429 {
 430     const int group[] = { 0, 1, 4, 5, 6, 7 };
 431     topManager_.initAtoms(9);
 432     topManager_.initUniformResidues(3);
 433
 434     gmx_ana_poscalc_t *pc1 = createCalculation(POS_RES, 0);
 435     gmx_ana_poscalc_t *pc2 = createCalculation(POS_RES, 0);
 436     gmx_ana_poscalc_t *pc3 = createCalculation(POS_RES, 0);
 437     setMaximumGroup(pc1, group);
 438     setMaximumGroup(pc2, group);
 439     setMaximumGroup(pc3, group);
 440     gmx_ana_pos_t *p1 = initPositions(pc1, "Positions");
 441     gmx_ana_pos_t *p2 = initPositions(pc2, "Positions");
 442     gmx_ana_pos_t *p3 = initPositions(pc3, "Positions");
 443     checkInitialized();
 444     {
 445         pcc_.initEvaluation();
 446         pcc_.initFrame();
 447         generateCoordinates();
 448         gmx::test::TestReferenceChecker frameCompound(
 449                 checker_.checkCompound("EvaluatedPositions", "Frame0"));
 450         updateAndCheck(pc1, p1, group, &frameCompound, "Positions");
 451         updateAndCheck(pc2, p2, group, &frameCompound, "Positions");
 452         updateAndCheck(pc3, p3, group, &frameCompound, "Positions");
 453     }
 454 }
 455
 456 TEST_F(PositionCalculationTest, HandlesOverlappingStaticCalculations)
 457 {
 458     const int group1[] = { 0, 1, 4, 5 };
 459     const int group2[] = { 4, 5, 7, 8 };
 460     topManager_.initAtoms(9);
 461     topManager_.initUniformResidues(3);
 462
 463     gmx_ana_poscalc_t *pc1 = createCalculation(POS_RES, 0);
 464     gmx_ana_poscalc_t *pc2 = createCalculation(POS_RES, 0);
 465     setMaximumGroup(pc1, group1);
 466     setMaximumGroup(pc2, group2);
 467     gmx_ana_pos_t *p1 = initPositions(pc1, "P1");
 468     gmx_ana_pos_t *p2 = initPositions(pc2, "P2");
 469     checkInitialized();
 470     {
 471         pcc_.initEvaluation();
 472         pcc_.initFrame();
 473         generateCoordinates();
 474         gmx::test::TestReferenceChecker frameCompound(
 475                 checker_.checkCompound("EvaluatedPositions", "Frame0"));
 476         updateAndCheck(pc1, p1, group1, &frameCompound, "P1");
 477         updateAndCheck(pc2, p2, group2, &frameCompound, "P2");
 478     }
 479 }
 480
 481 // TODO: Check for handling of more multiple calculation cases
 482
 483 } // namespace