src/gromacs/trajectoryanalysis/tests/topologyinformation.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 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 for TopologyInformation
  38  *
  39  * \author Mark Abraham <mark.j.abraham@gmail.com>
  40  * \ingroup module_trajectoryanalysis
  41  */
  42 #include "gmxpre.h"
  43
  44 #include "gromacs/trajectoryanalysis/topologyinformation.h"
  45
  46 #include <gtest/gtest.h>
  47
  48 #include "gromacs/gmxpreprocess/grompp.h"
  49 #include "gromacs/math/vectypes.h"
  50 #include "gromacs/topology/topology.h"
  51 #include "gromacs/utility/arrayref.h"
  52 #include "gromacs/utility/exceptions.h"
  53 #include "gromacs/utility/textwriter.h"
  54
  55 #include "testutils/cmdlinetest.h"
  56 #include "testutils/testfilemanager.h"
  57
  58 #include "moduletest.h"
  59
  60 namespace gmx
  61 {
  62 namespace test
  63 {
  64 namespace
  65 {
  66
  67 TEST(TopologyInformation, CantWorkWithoutReadingAFile)
  68 {
  69     TopologyInformation topInfo;
  70     EXPECT_FALSE(topInfo.hasTopology());
  71     EXPECT_FALSE(topInfo.hasFullTopology());
  72     EXPECT_EQ(nullptr, topInfo.mtop());
  73     EXPECT_EQ(nullptr, topInfo.expandedTopology());
  74     auto atoms1 = topInfo.copyAtoms();
  75     EXPECT_TRUE(atoms1);
  76     auto atoms2 = topInfo.copyAtoms();
  77     ASSERT_TRUE(atoms2);
  78     EXPECT_NE(atoms1.get(), atoms2.get());
  79     EXPECT_EQ(0, atoms1->nr);
  80     EXPECT_EQ(-1, topInfo.ePBC());
  81     EXPECT_THROW(topInfo.x().size(), gmx::APIError);
  82     EXPECT_THROW(topInfo.v().size(), gmx::APIError);
  83     matrix box{ { -2 } };
  84     topInfo.getBox(box);
  85     EXPECT_EQ(0, box[XX][XX]);
  86     EXPECT_EQ(0, box[XX][YY]);
  87     EXPECT_EQ(0, box[XX][ZZ]);
  88     EXPECT_EQ(0, box[YY][XX]);
  89     EXPECT_EQ(0, box[YY][YY]);
  90     EXPECT_EQ(0, box[YY][ZZ]);
  91     EXPECT_EQ(0, box[ZZ][XX]);
  92     EXPECT_EQ(0, box[ZZ][YY]);
  93     EXPECT_EQ(0, box[ZZ][ZZ]);
  94     EXPECT_FALSE(topInfo.name());
  95 }
  96
  97 //! Common test code to reduce duplication
  98 void runCommonTests(const TopologyInformation& topInfo, const int numAtoms)
  99 {
 100     EXPECT_TRUE(topInfo.hasTopology());
 101     ASSERT_TRUE(topInfo.mtop());
 102     EXPECT_EQ(numAtoms, topInfo.mtop()->natoms);
 103     // TODO Dump mtop to refdata when that is possible
 104     ASSERT_TRUE(topInfo.expandedTopology());
 105     auto atoms1 = topInfo.copyAtoms();
 106     EXPECT_TRUE(atoms1);
 107     auto atoms2 = topInfo.copyAtoms();
 108     EXPECT_TRUE(atoms2);
 109     // Must be different pointer to a deep copy.
 110     EXPECT_NE(atoms1.get(), atoms2.get());
 111     auto atoms = topInfo.atoms();
 112     // Must be a pointer to a different instance.
 113     EXPECT_NE(atoms1.get(), atoms);
 114     EXPECT_NE(atoms2.get(), atoms);
 115     EXPECT_EQ(numAtoms, topInfo.x().size());
 116     EXPECT_EQ(numAtoms, topInfo.v().size());
 117     matrix box{ { -2 } };
 118     topInfo.getBox(box);
 119     EXPECT_FLOAT_EQ(5.9062, box[XX][XX]);
 120     EXPECT_FLOAT_EQ(0, box[XX][YY]);
 121     EXPECT_FLOAT_EQ(0, box[XX][ZZ]);
 122     EXPECT_FLOAT_EQ(0, box[YY][XX]);
 123     EXPECT_FLOAT_EQ(6.8451, box[YY][YY]);
 124     EXPECT_FLOAT_EQ(0, box[YY][ZZ]);
 125     EXPECT_FLOAT_EQ(0, box[ZZ][XX]);
 126     EXPECT_FLOAT_EQ(0, box[ZZ][YY]);
 127     EXPECT_FLOAT_EQ(3.0517, box[ZZ][ZZ]);
 128     EXPECT_STREQ("First 10 residues from 1AKI", topInfo.name());
 129 }
 130
 131 TEST(TopologyInformation, WorksWithGroFile)
 132 {
 133     const int           numAtoms = 156;
 134     TopologyInformation topInfo;
 135     topInfo.fillFromInputFile(TestFileManager::getInputFilePath("lysozyme.gro"));
 136     EXPECT_FALSE(topInfo.hasFullTopology());
 137     runCommonTests(topInfo, numAtoms);
 138     EXPECT_EQ(-1, topInfo.ePBC());
 139
 140     // Check the per-atom data
 141     auto atoms = topInfo.copyAtoms();
 142     ASSERT_EQ(numAtoms, atoms->nr);
 143     EXPECT_TRUE(atoms->haveMass);
 144     // TODO atommass.dat assumes united atom CA, which is probably not expected behaviour
 145     EXPECT_FLOAT_EQ(13.019, atoms->atom[26].m);
 146     EXPECT_FALSE(atoms->haveCharge);
 147     EXPECT_FALSE(atoms->haveType);
 148     EXPECT_EQ(0, atoms->atom[26].type);
 149     EXPECT_EQ(0, atoms->atom[26].atomnumber);
 150     EXPECT_EQ(1, atoms->atom[26].resind);
 151     // gro files don't have the information that pdb files might
 152     EXPECT_FALSE(atoms->havePdbInfo);
 153     EXPECT_FALSE(atoms->pdbinfo);
 154     EXPECT_EQ(10, atoms->nres);
 155     ASSERT_TRUE(atoms->resinfo);
 156     ASSERT_TRUE(atoms->resinfo[4].name);
 157     EXPECT_STREQ("ARG", *atoms->resinfo[4].name);
 158     EXPECT_EQ(5, atoms->resinfo[4].nr);
 159     EXPECT_EQ(0, atoms->resinfo[4].chainnum);
 160     EXPECT_EQ(' ', atoms->resinfo[4].chainid);
 161 }
 162
 163 TEST(TopologyInformation, WorksWithPdbFile)
 164 {
 165     const int           numAtoms = 156;
 166     TopologyInformation topInfo;
 167     topInfo.fillFromInputFile(TestFileManager::getInputFilePath("lysozyme.pdb"));
 168     EXPECT_FALSE(topInfo.hasFullTopology());
 169     runCommonTests(topInfo, numAtoms);
 170     // TODO why does this differ from .gro?
 171     EXPECT_EQ(0, topInfo.ePBC());
 172
 173     // Check the per-atom data
 174     auto atoms = topInfo.copyAtoms();
 175     ASSERT_EQ(numAtoms, atoms->nr);
 176     EXPECT_TRUE(atoms->haveMass);
 177     // TODO atommass.dat assumes united atom CA, which is probably not expected behaviour
 178     EXPECT_FLOAT_EQ(13.019, atoms->atom[26].m);
 179     EXPECT_FALSE(atoms->haveCharge);
 180     EXPECT_FALSE(atoms->haveType);
 181     EXPECT_EQ(0, atoms->atom[26].type);
 182     EXPECT_EQ(0, atoms->atom[26].atomnumber);
 183     EXPECT_EQ(1, atoms->atom[26].resind);
 184     // pdb files can carry more information than gro
 185     EXPECT_TRUE(atoms->havePdbInfo);
 186     ASSERT_TRUE(atoms->pdbinfo);
 187     EXPECT_EQ(10, atoms->nres);
 188     ASSERT_TRUE(atoms->resinfo);
 189     ASSERT_TRUE(atoms->resinfo[4].name);
 190     EXPECT_STREQ("ARG", *atoms->resinfo[4].name);
 191     EXPECT_EQ(5, atoms->resinfo[4].nr);
 192     EXPECT_EQ(0, atoms->resinfo[4].chainnum);
 193     EXPECT_EQ('B', atoms->resinfo[4].chainid);
 194 }
 195
 196 TEST(TopologyInformation, WorksWithTprFromPdbFile)
 197 {
 198     TestFileManager fileManager;
 199
 200     // Make the tpr file to use
 201     std::string       name             = "lysozyme";
 202     const std::string mdpInputFileName = fileManager.getTemporaryFilePath(name + ".mdp");
 203     // Ensure the seeds have a value so that the resulting .tpr dump
 204     // is reproducible.
 205     TextWriter::writeFileFromString(mdpInputFileName, "");
 206     std::string tprName = fileManager.getTemporaryFilePath(name + ".tpr");
 207     {
 208         CommandLine caller;
 209         caller.append("grompp");
 210         caller.addOption("-f", mdpInputFileName);
 211         caller.addOption("-p", TestFileManager::getInputFilePath(name));
 212         caller.addOption("-c", TestFileManager::getInputFilePath(name + ".pdb"));
 213         caller.addOption("-o", tprName);
 214         ASSERT_EQ(0, gmx_grompp(caller.argc(), caller.argv()));
 215     }
 216
 217     const int           numAtoms = 156;
 218     TopologyInformation topInfo;
 219     topInfo.fillFromInputFile(tprName);
 220     EXPECT_TRUE(topInfo.hasFullTopology());
 221     runCommonTests(topInfo, numAtoms);
 222     // TODO why does this differ from .gro?
 223     EXPECT_EQ(0, topInfo.ePBC());
 224
 225     // Check the per-atom data
 226     auto atoms = topInfo.copyAtoms();
 227     ASSERT_EQ(numAtoms, atoms->nr);
 228     EXPECT_TRUE(atoms->haveMass);
 229     EXPECT_FLOAT_EQ(12.011, atoms->atom[26].m);
 230     EXPECT_TRUE(atoms->haveCharge);
 231     EXPECT_TRUE(atoms->haveType);
 232     EXPECT_EQ(2, atoms->atom[26].type);
 233     EXPECT_EQ(6, atoms->atom[26].atomnumber);
 234     EXPECT_EQ(1, atoms->atom[26].resind);
 235     // tpr files also don't carry pdb information
 236     EXPECT_FALSE(atoms->havePdbInfo);
 237     EXPECT_FALSE(atoms->pdbinfo);
 238     EXPECT_EQ(10, atoms->nres);
 239     ASSERT_TRUE(atoms->resinfo);
 240     ASSERT_TRUE(atoms->resinfo[4].name);
 241     EXPECT_STREQ("ARG", *atoms->resinfo[4].name);
 242     EXPECT_EQ(5, atoms->resinfo[4].nr);
 243     EXPECT_EQ(0, atoms->resinfo[4].chainnum);
 244     // In particular, chain ID does not get recorded in the .tpr file
 245     EXPECT_EQ(0, atoms->resinfo[4].chainid);
 246 }
 247
 248 } // namespace
 249 } // namespace test
 250 } // namespace gmx