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