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