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