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.
   5  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   6  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   7  * and including many others, as listed in the AUTHORS file in the
   8  * top-level source directory and at http://www.gromacs.org.
   9  *
  10  * GROMACS is free software; you can redistribute it and/or
  11  * modify it under the terms of the GNU Lesser General Public License
  12  * as published by the Free Software Foundation; either version 2.1
  13  * of the License, or (at your option) any later version.
  14  *
  15  * GROMACS is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18  * Lesser General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU Lesser General Public
  21  * License along with GROMACS; if not, see
  22  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  23  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  24  *
  25  * If you want to redistribute modifications to GROMACS, please
  26  * consider that scientific software is very special. Version
  27  * control is crucial - bugs must be traceable. We will be happy to
  28  * consider code for inclusion in the official distribution, but
  29  * derived work must not be called official GROMACS. Details are found
  30  * in the README & COPYING files - if they are missing, get the
  31  * official version at http://www.gromacs.org.
  32  *
  33  * To help us fund GROMACS development, we humbly ask that you cite
  34  * the research papers on the package. Check out http://www.gromacs.org.
  35  */
  36 /*! \internal \file
  37  * \brief
  38  * Tests the position mapping engine.
  39  *
  40  * \author Teemu Murtola <teemu.murtola@gmail.com>
  41  * \ingroup module_selection
  42  */
  43 #include "gmxpre.h"
  44
  45 #include "gromacs/selection/poscalc.h"
  46
  47 #include <memory>
  48 #include <vector>
  49
  50 #include <gtest/gtest.h>
  51
  52 #include "gromacs/math/vec.h"
  53 #include "gromacs/selection/indexutil.h"
  54 #include "gromacs/selection/position.h"
  55 #include "gromacs/topology/topology.h"
  56 #include "gromacs/trajectory/trajectoryframe.h"
  57 #include "gromacs/utility/arrayref.h"
  58 #include "gromacs/utility/smalloc.h"
  59
  60 #include "testutils/refdata.h"
  61
  62 #include "toputils.h"
  63
  64 namespace
  65 {
  66
  67 /********************************************************************
  68  * PositionCalculationTest
  69  */
  70
  71 class PositionCalculationTest : public ::testing::Test
  72 {
  73 public:
  74     PositionCalculationTest();
  75     ~PositionCalculationTest() override;
  76
  77     void generateCoordinates();
  78
  79     gmx_ana_poscalc_t* createCalculation(e_poscalc_t type, int flags);
  80     void           setMaximumGroup(gmx_ana_poscalc_t* pc, const gmx::ArrayRef<const int>& atoms);
  81     gmx_ana_pos_t* initPositions(gmx_ana_poscalc_t* pc, const char* name);
  82
  83     void checkInitialized();
  84     void updateAndCheck(gmx_ana_poscalc_t*               pc,
  85                         gmx_ana_pos_t*                   p,
  86                         const gmx::ArrayRef<const int>&  atoms,
  87                         gmx::test::TestReferenceChecker* checker,
  88                         const char*                      name);
  89
  90     void testSingleStatic(e_poscalc_t                     type,
  91                           int                             flags,
  92                           bool                            bExpectTop,
  93                           const gmx::ArrayRef<const int>& atoms,
  94                           const gmx::ArrayRef<const int>& index = {});
  95     void testSingleDynamic(e_poscalc_t                     type,
  96                            int                             flags,
  97                            bool                            bExpectTop,
  98                            const gmx::ArrayRef<const int>& initAtoms,
  99                            const gmx::ArrayRef<const int>& evalAtoms,
 100                            const gmx::ArrayRef<const int>& index = {});
 101
 102     gmx::test::TestReferenceData       data_;
 103     gmx::test::TestReferenceChecker    checker_;
 104     gmx::test::TopologyManager         topManager_;
 105     gmx::PositionCalculationCollection pcc_;
 106
 107 private:
 108     typedef std::unique_ptr<gmx_ana_pos_t> PositionPointer;
 109
 110     struct PositionTest
 111     {
 112         PositionTest(PositionPointer pos, gmx_ana_poscalc_t* pc, const char* name) :
 113             pos(std::move(pos)),
 114             pc(pc),
 115             name(name)
 116         {
 117         }
 118
 119         PositionPointer    pos;
 120         gmx_ana_poscalc_t* pc;
 121         const char*        name;
 122     };
 123
 124     typedef std::vector<PositionTest> PositionTestList;
 125
 126     void setTopologyIfRequired();
 127     void checkPositions(gmx::test::TestReferenceChecker* checker, const char* name, gmx_ana_pos_t* p, bool bCoordinates);
 128
 129     std::vector<gmx_ana_poscalc_t*> pcList_;
 130     PositionTestList                posList_;
 131     bool                            bTopSet_;
 132 };
 133
 134 PositionCalculationTest::PositionCalculationTest() : checker_(data_.rootChecker()), bTopSet_(false)
 135 {
 136     topManager_.requestFrame();
 137 }
 138
 139 PositionCalculationTest::~PositionCalculationTest()
 140 {
 141     std::vector<gmx_ana_poscalc_t*>::reverse_iterator pci;
 142     for (pci = pcList_.rbegin(); pci != pcList_.rend(); ++pci)
 143     {
 144         gmx_ana_poscalc_free(*pci);
 145     }
 146 }
 147
 148 void PositionCalculationTest::generateCoordinates()
 149 {
 150     t_atoms&    atoms = topManager_.atoms();
 151     t_trxframe* frame = topManager_.frame();
 152     for (int i = 0; i < atoms.nr; ++i)
 153     {
 154         frame->x[i][XX] = i;
 155         frame->x[i][YY] = atoms.atom[i].resind;
 156         frame->x[i][ZZ] = 0.0;
 157         if (frame->bV)
 158         {
 159             copy_rvec(frame->x[i], frame->v[i]);
 160             frame->v[i][ZZ] = 1.0;
 161         }
 162         if (frame->bF)
 163         {
 164             copy_rvec(frame->x[i], frame->f[i]);
 165             frame->f[i][ZZ] = -1.0;
 166         }
 167     }
 168 }
 169
 170 gmx_ana_poscalc_t* PositionCalculationTest::createCalculation(e_poscalc_t type, int flags)
 171 {
 172     pcList_.reserve(pcList_.size() + 1);
 173     pcList_.push_back(pcc_.createCalculation(type, flags));
 174     return pcList_.back();
 175 }
 176
 177 void PositionCalculationTest::setMaximumGroup(gmx_ana_poscalc_t* pc, const gmx::ArrayRef<const int>& atoms)
 178 {
 179     setTopologyIfRequired();
 180     gmx_ana_index_t g;
 181     g.isize = atoms.size();
 182     g.index = const_cast<int*>(atoms.data());
 183     gmx_ana_poscalc_set_maxindex(pc, &g);
 184 }
 185
 186 gmx_ana_pos_t* PositionCalculationTest::initPositions(gmx_ana_poscalc_t* pc, const char* name)
 187 {
 188     posList_.reserve(posList_.size() + 1);
 189     PositionPointer p(new gmx_ana_pos_t());
 190     gmx_ana_pos_t*  result = p.get();
 191     posList_.emplace_back(std::move(p), pc, name);
 192     gmx_ana_poscalc_init_pos(pc, result);
 193     return result;
 194 }
 195
 196 void PositionCalculationTest::checkInitialized()
 197 {
 198     gmx::test::TestReferenceChecker compound(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(gmx_ana_poscalc_t*               pc,
 207                                              gmx_ana_pos_t*                   p,
 208                                              const gmx::ArrayRef<const int>&  atoms,
 209                                              gmx::test::TestReferenceChecker* checker,
 210                                              const char*                      name)
 211 {
 212     gmx_ana_index_t g;
 213     g.isize = atoms.size();
 214     g.index = const_cast<int*>(atoms.data());
 215     gmx_ana_poscalc_update(pc, p, &g, topManager_.frame(), nullptr);
 216     checkPositions(checker, name, p, true);
 217 }
 218
 219 void PositionCalculationTest::testSingleStatic(e_poscalc_t                     type,
 220                                                int                             flags,
 221                                                bool                            bExpectTop,
 222                                                const gmx::ArrayRef<const int>& atoms,
 223                                                const gmx::ArrayRef<const int>& index)
 224 {
 225     t_trxframe* frame = topManager_.frame();
 226     if (frame->bV)
 227     {
 228         flags |= POS_VELOCITIES;
 229     }
 230     if (frame->bF)
 231     {
 232         flags |= POS_FORCES;
 233     }
 234     gmx_ana_poscalc_t* pc               = createCalculation(type, flags);
 235     const bool         requiresTopology = gmx_ana_poscalc_required_topology_info(pc)
 236                                   != gmx::PositionCalculationCollection::RequiredTopologyInfo::None;
 237     EXPECT_EQ(bExpectTop, requiresTopology);
 238     setMaximumGroup(pc, atoms);
 239     gmx_ana_pos_t* p = initPositions(pc, nullptr);
 240     checkInitialized();
 241     {
 242         generateCoordinates();
 243         if (!index.empty())
 244         {
 245             topManager_.initFrameIndices(index);
 246         }
 247         pcc_.initEvaluation();
 248         pcc_.initFrame(frame);
 249         gmx::test::TestReferenceChecker frameCompound(
 250                 checker_.checkCompound("EvaluatedPositions", "Frame0"));
 251         updateAndCheck(pc, p, atoms, &frameCompound, nullptr);
 252     }
 253 }
 254
 255 void PositionCalculationTest::testSingleDynamic(e_poscalc_t                     type,
 256                                                 int                             flags,
 257                                                 bool                            bExpectTop,
 258                                                 const gmx::ArrayRef<const int>& initAtoms,
 259                                                 const gmx::ArrayRef<const int>& evalAtoms,
 260                                                 const gmx::ArrayRef<const int>& index)
 261 {
 262     gmx_ana_poscalc_t* pc               = createCalculation(type, flags | POS_DYNAMIC);
 263     const bool         requiresTopology = gmx_ana_poscalc_required_topology_info(pc)
 264                                   != gmx::PositionCalculationCollection::RequiredTopologyInfo::None;
 265     EXPECT_EQ(bExpectTop, requiresTopology);
 266     setMaximumGroup(pc, initAtoms);
 267     gmx_ana_pos_t* p = initPositions(pc, nullptr);
 268     checkInitialized();
 269     {
 270         generateCoordinates();
 271         if (!index.empty())
 272         {
 273             topManager_.initFrameIndices(index);
 274         }
 275         pcc_.initEvaluation();
 276         pcc_.initFrame(topManager_.frame());
 277         gmx::test::TestReferenceChecker frameCompound(
 278                 checker_.checkCompound("EvaluatedPositions", "Frame0"));
 279         updateAndCheck(pc, p, evalAtoms, &frameCompound, nullptr);
 280     }
 281 }
 282
 283 void PositionCalculationTest::setTopologyIfRequired()
 284 {
 285     if (bTopSet_)
 286     {
 287         return;
 288     }
 289     std::vector<gmx_ana_poscalc_t*>::const_iterator pci;
 290     for (pci = pcList_.begin(); pci != pcList_.end(); ++pci)
 291     {
 292         const bool requiresTopology = gmx_ana_poscalc_required_topology_info(*pci)
 293                                       != gmx::PositionCalculationCollection::RequiredTopologyInfo::None;
 294         if (requiresTopology)
 295         {
 296             bTopSet_ = true;
 297             pcc_.setTopology(topManager_.topology());
 298             return;
 299         }
 300     }
 301 }
 302
 303 void PositionCalculationTest::checkPositions(gmx::test::TestReferenceChecker* checker,
 304                                              const char*                      name,
 305                                              gmx_ana_pos_t*                   p,
 306                                              bool                             bCoordinates)
 307 {
 308     gmx::test::TestReferenceChecker compound(checker->checkCompound("Positions", name));
 309     compound.checkInteger(p->count(), "Count");
 310     const char* type = "???";
 311     switch (p->m.type)
 312     {
 313         case INDEX_UNKNOWN: type = "unknown"; break;
 314         case INDEX_ATOM: type = "atoms"; break;
 315         case INDEX_RES: type = "residues"; break;
 316         case INDEX_MOL: type = "molecules"; break;
 317         case INDEX_ALL: type = "single"; break;
 318     }
 319     compound.checkString(type, "Type");
 320     compound.checkSequenceArray(p->count() + 1, p->m.mapb.index, "Block");
 321     for (int i = 0; i < p->count(); ++i)
 322     {
 323         gmx::test::TestReferenceChecker posCompound(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]], "Atoms");
 326         posCompound.checkInteger(p->m.refid[i], "RefId");
 327         if (bCoordinates)
 328         {
 329             posCompound.checkVector(p->x[i], "Coordinates");
 330         }
 331         if (bCoordinates && p->v != nullptr)
 332         {
 333             posCompound.checkVector(p->v[i], "Velocity");
 334         }
 335         if (bCoordinates && p->f != nullptr)
 336         {
 337             posCompound.checkVector(p->f[i], "Force");
 338         }
 339         int originalIdIndex = (p->m.refid[i] != -1 ? p->m.refid[i] : i);
 340         EXPECT_EQ(p->m.orgid[originalIdIndex], p->m.mapid[i]);
 341     }
 342 }
 343
 344 /********************************************************************
 345  * Actual tests
 346  */
 347
 348 TEST_F(PositionCalculationTest, ComputesAtomPositions)
 349 {
 350     const int group[] = { 0, 1, 2, 3 };
 351     topManager_.requestVelocities();
 352     topManager_.requestForces();
 353     topManager_.initAtoms(4);
 354     testSingleStatic(POS_ATOM, 0, false, group);
 355 }
 356
 357 TEST_F(PositionCalculationTest, ComputesResidueCOGPositions)
 358 {
 359     const int group[] = { 0, 1, 2, 3, 4, 8 };
 360     topManager_.requestVelocities();
 361     topManager_.requestForces();
 362     topManager_.initAtoms(9);
 363     topManager_.initUniformResidues(3);
 364     testSingleStatic(POS_RES, 0, true, group);
 365 }
 366
 367 TEST_F(PositionCalculationTest, ComputesResidueCOMPositions)
 368 {
 369     const int group[] = { 0, 1, 2, 3, 4, 8 };
 370     topManager_.requestVelocities();
 371     topManager_.requestForces();
 372     topManager_.initAtoms(9);
 373     topManager_.initUniformResidues(3);
 374     testSingleStatic(POS_RES, POS_MASS, true, group);
 375 }
 376
 377 TEST_F(PositionCalculationTest, ComputesGroupCOGPositions)
 378 {
 379     const int group[] = { 0, 1, 2, 3, 4, 5, 6, 7, 8 };
 380     topManager_.requestVelocities();
 381     topManager_.requestForces();
 382     topManager_.initAtoms(9);
 383     // Topology (masses) is requires for computing the force
 384     testSingleStatic(POS_ALL, 0, true, group);
 385 }
 386
 387 TEST_F(PositionCalculationTest, ComputesGroupCOMPositions)
 388 {
 389     const int group[] = { 0, 1, 2, 3, 4, 5, 6, 7, 8 };
 390     topManager_.requestVelocities();
 391     topManager_.requestForces();
 392     topManager_.initAtoms(9);
 393     testSingleStatic(POS_ALL, POS_MASS, true, group);
 394 }
 395
 396 TEST_F(PositionCalculationTest, ComputesPositionsWithCompleteWhole)
 397 {
 398     const int group[] = { 0, 1, 2, 3, 4, 8 };
 399     topManager_.initAtoms(9);
 400     topManager_.initUniformResidues(3);
 401     testSingleStatic(POS_RES, POS_COMPLWHOLE, true, group);
 402 }
 403
 404 TEST_F(PositionCalculationTest, ComputesPositionsWithCompleteMax)
 405 {
 406     const int maxGroup[]  = { 0, 1, 4, 5, 6, 8 };
 407     const int evalGroup[] = { 0, 1, 5, 6 };
 408     topManager_.initAtoms(9);
 409     topManager_.initUniformResidues(3);
 410     testSingleDynamic(POS_RES, POS_COMPLMAX, true, maxGroup, evalGroup);
 411 }
 412
 413 TEST_F(PositionCalculationTest, ComputesPositionMask)
 414 {
 415     const int maxGroup[]  = { 0, 1, 2, 3, 4, 5 };
 416     const int evalGroup[] = { 1, 2, 4 };
 417     topManager_.initAtoms(6);
 418     testSingleDynamic(POS_ATOM, POS_MASKONLY, false, maxGroup, evalGroup);
 419 }
 420
 421 // TODO: Check for POS_ALL_PBC
 422
 423 TEST_F(PositionCalculationTest, HandlesFramesWithLessAtoms)
 424 {
 425     const int group[] = { 2, 3, 5, 6 };
 426     topManager_.initAtoms(10);
 427     testSingleStatic(POS_ATOM, 0, false, group, group);
 428 }
 429
 430 TEST_F(PositionCalculationTest, HandlesFramesWithLessAtoms2)
 431 {
 432     const int group[] = { 2, 3, 6, 7 };
 433     const int index[] = { 1, 2, 3, 4, 6, 7 };
 434     topManager_.initAtoms(10);
 435     testSingleStatic(POS_ATOM, 0, false, group, index);
 436 }
 437
 438 TEST_F(PositionCalculationTest, HandlesIdenticalStaticCalculations)
 439 {
 440     const int group[] = { 0, 1, 4, 5, 6, 7 };
 441     topManager_.initAtoms(9);
 442     topManager_.initUniformResidues(3);
 443
 444     gmx_ana_poscalc_t* pc1 = createCalculation(POS_RES, 0);
 445     gmx_ana_poscalc_t* pc2 = createCalculation(POS_RES, 0);
 446     gmx_ana_poscalc_t* pc3 = createCalculation(POS_RES, 0);
 447     setMaximumGroup(pc1, group);
 448     setMaximumGroup(pc2, group);
 449     setMaximumGroup(pc3, group);
 450     gmx_ana_pos_t* p1 = initPositions(pc1, "Positions");
 451     gmx_ana_pos_t* p2 = initPositions(pc2, "Positions");
 452     gmx_ana_pos_t* p3 = initPositions(pc3, "Positions");
 453     checkInitialized();
 454     {
 455         generateCoordinates();
 456         pcc_.initEvaluation();
 457         pcc_.initFrame(topManager_.frame());
 458         gmx::test::TestReferenceChecker frameCompound(
 459                 checker_.checkCompound("EvaluatedPositions", "Frame0"));
 460         updateAndCheck(pc1, p1, group, &frameCompound, "Positions");
 461         updateAndCheck(pc2, p2, group, &frameCompound, "Positions");
 462         updateAndCheck(pc3, p3, group, &frameCompound, "Positions");
 463     }
 464 }
 465
 466 TEST_F(PositionCalculationTest, HandlesOverlappingStaticCalculations)
 467 {
 468     const int group1[] = { 0, 1, 4, 5 };
 469     const int group2[] = { 4, 5, 7, 8 };
 470     topManager_.initAtoms(9);
 471     topManager_.initUniformResidues(3);
 472
 473     gmx_ana_poscalc_t* pc1 = createCalculation(POS_RES, 0);
 474     gmx_ana_poscalc_t* pc2 = createCalculation(POS_RES, 0);
 475     setMaximumGroup(pc1, group1);
 476     setMaximumGroup(pc2, group2);
 477     gmx_ana_pos_t* p1 = initPositions(pc1, "P1");
 478     gmx_ana_pos_t* p2 = initPositions(pc2, "P2");
 479     checkInitialized();
 480     {
 481         generateCoordinates();
 482         pcc_.initEvaluation();
 483         pcc_.initFrame(topManager_.frame());
 484         gmx::test::TestReferenceChecker frameCompound(
 485                 checker_.checkCompound("EvaluatedPositions", "Frame0"));
 486         updateAndCheck(pc1, p1, group1, &frameCompound, "P1");
 487         updateAndCheck(pc2, p2, group2, &frameCompound, "P2");
 488     }
 489 }
 490
 491 // TODO: Check for handling of more multiple calculation cases
 492
 493 } // namespace