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