src/gromacs/trajectoryanalysis/tests/topologyinformation.cpp

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