src/gromacs/selection/tests/poscalc.cpp

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