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