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