src/gromacs/applied_forces/qmmm/tests/qmmmtopologypreprocessor.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2021, 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 QMMMInputGenerator class for QMMM MDModule
  38  *
  39  * \author Dmitry Morozov <dmitry.morozov@jyu.fi>
  40  * \ingroup module_applied_forces
  41  */
  42 #include "gmxpre.h"
  43
  44 #include "gromacs/applied_forces/qmmm/qmmmtopologypreprocessor.h"
  45
  46 #include <vector>
  47
  48 #include <gtest/gtest.h>
  49
  50 #include "gromacs/fileio/confio.h"
  51 #include "gromacs/gmxpreprocess/grompp.h"
  52 #include "gromacs/math/vec.h"
  53 #include "gromacs/pbcutil/pbc.h"
  54 #include "gromacs/topology/mtop_lookup.h"
  55 #include "gromacs/topology/mtop_util.h"
  56 #include "gromacs/topology/topology.h"
  57 #include "gromacs/utility/arrayref.h"
  58 #include "gromacs/utility/path.h"
  59 #include "gromacs/utility/textwriter.h"
  60 #include "testutils/cmdlinetest.h"
  61 #include "testutils/refdata.h"
  62 #include "testutils/testasserts.h"
  63 #include "testutils/testfilemanager.h"
  64
  65 namespace gmx
  66 {
  67
  68 class QMMMTopologyPreprocessorTest : public ::testing::Test
  69 {
  70 public:
  71     /*! \brief Generates tpr file from *.top and *.gro existing in the simulation database directory
  72      * and loads gmx_mtop_t from it
  73      */
  74     void makeMtopFromFile(const std::string& fileName, const std::string& mdpContent)
  75     {
  76         const std::string simData = gmx::test::TestFileManager::getTestSimulationDatabaseDirectory();
  77
  78         // Generate empty mdp file
  79         const std::string mdpInputFileName = fileManager_.getTemporaryFilePath(fileName + ".mdp");
  80         gmx::TextWriter::writeFileFromString(mdpInputFileName, mdpContent);
  81
  82         // Generate tpr file
  83         const std::string tprName = fileManager_.getTemporaryFilePath(fileName + ".tpr");
  84         {
  85             gmx::test::CommandLine caller;
  86             caller.append("grompp");
  87             caller.addOption("-f", mdpInputFileName);
  88             caller.addOption("-p", gmx::Path::join(simData, fileName + ".top"));
  89             caller.addOption("-c", gmx::Path::join(simData, fileName + ".gro"));
  90             caller.addOption("-o", tprName);
  91             ASSERT_EQ(0, gmx_grompp(caller.argc(), caller.argv()));
  92         }
  93
  94         // Load topology
  95         bool    fullTopology;
  96         PbcType pbcType;
  97         matrix  box;
  98         readConfAndTopology(tprName.c_str(), &fullTopology, &mtop_, &pbcType, nullptr, nullptr, box);
  99     }
 100
 101 protected:
 102     gmx::test::TestFileManager fileManager_;
 103     std::vector<index>         qmIndices_;
 104     gmx_mtop_t                 mtop_;
 105 };
 106
 107 class QMMMTopologyPreprocessorChecker : public gmx::test::TestReferenceChecker
 108 {
 109 public:
 110     QMMMTopologyPreprocessorChecker(const gmx::test::TestReferenceChecker& rootChecker) :
 111         gmx::test::TestReferenceChecker(rootChecker)
 112     {
 113         // Tolerance of all charges and vectors should be 1E-3
 114         setDefaultTolerance(gmx::test::absoluteTolerance(0.001));
 115     }
 116
 117     void checkAll(const QMMMTopologyInfo& info, const gmx_mtop_t& mtop)
 118     {
 119         checkInteger(info.numMMAtoms, "Number of MM atoms");
 120         checkInteger(info.numQMAtoms, "Number of QM atoms");
 121         checkReal(info.remainingMMCharge, "MM charge");
 122         checkReal(info.totalClassicalChargeOfQMAtoms, "QM charge");
 123         checkInteger(info.numExclusionsMade, "Exclusions Made");
 124         checkInteger(info.numBondsRemoved, "Removed bonds");
 125         checkInteger(info.numAnglesRemoved, "Removed angles");
 126         checkInteger(info.numDihedralsRemoved, "Removed dihedrals");
 127         checkInteger(info.numSettleRemoved, "Removed settles");
 128         checkInteger(info.numConnBondsAdded, "Generated CONNBONDS");
 129         checkInteger(info.numVirtualSitesModified, "Removed vsites");
 130         checkInteger(info.numConstrainedBondsInQMSubsystem, "QM Constraints Found");
 131         checkInteger(info.numLinkBonds, "Link-atom sites");
 132         checkInteger(mtop.molblock.size(), "Molblocks in topology");
 133     }
 134 };
 135
 136 TEST_F(QMMMTopologyPreprocessorTest, CanConstruct)
 137 {
 138     qmIndices_ = { 0, 1, 2 };
 139     EXPECT_NO_THROW(QMMMTopologyPreprocessor topPrep(qmIndices_));
 140 }
 141
 142 TEST_F(QMMMTopologyPreprocessorTest, FourWatersFirstQMNoLink)
 143 {
 144     // Reference input 4x SPCE waters from database 4waters.top (first one QM) and no Link atoms
 145     makeMtopFromFile("4water", "");
 146     qmIndices_ = { 0, 1, 2 };
 147
 148     QMMMTopologyPreprocessor topPrep(qmIndices_);
 149     topPrep.preprocess(&mtop_);
 150
 151     // Get data about changes and check it
 152     QMMMTopologyInfo info = topPrep.topInfo();
 153
 154     gmx::test::TestReferenceData    data;
 155     QMMMTopologyPreprocessorChecker checker(data.rootChecker());
 156     checker.checkAll(info, mtop_);
 157 }
 158
 159 TEST_F(QMMMTopologyPreprocessorTest, FourWatersSeondAndForthQMNoLink)
 160 {
 161     // Reference input 4x SPCE waters from database 4waters.top (second and forth are QM) and no Link atoms
 162     makeMtopFromFile("4water", "");
 163     qmIndices_ = { 3, 4, 5, 9, 10, 11 };
 164
 165     QMMMTopologyPreprocessor topPrep(qmIndices_);
 166     topPrep.preprocess(&mtop_);
 167
 168     // Get data about changes and check it
 169     QMMMTopologyInfo info = topPrep.topInfo();
 170
 171     gmx::test::TestReferenceData    data;
 172     QMMMTopologyPreprocessorChecker checker(data.rootChecker());
 173     checker.checkAll(info, mtop_);
 174 }
 175
 176 TEST_F(QMMMTopologyPreprocessorTest, FourWatersFirstQMWithLink)
 177 {
 178     // Reference input 4x SPCE waters from database 4waters.top (first one QM) with Link atom
 179     makeMtopFromFile("4water", "");
 180     qmIndices_ = { 0, 1 };
 181
 182     QMMMTopologyPreprocessor topPrep(qmIndices_);
 183     topPrep.preprocess(&mtop_);
 184
 185     // Get data about changes and check it
 186     QMMMTopologyInfo info = topPrep.topInfo();
 187
 188     gmx::test::TestReferenceData    data;
 189     QMMMTopologyPreprocessorChecker checker(data.rootChecker());
 190     checker.checkAll(info, mtop_);
 191 }
 192
 193 TEST_F(QMMMTopologyPreprocessorTest, AlanineDipeptideWithLinksNoConstraints)
 194 {
 195     // Reference input alanine_vacuo.top
 196     makeMtopFromFile("alanine_vacuo", "");
 197     qmIndices_ = { 8, 9, 10, 11, 12, 13 };
 198
 199     QMMMTopologyPreprocessor topPrep(qmIndices_);
 200     topPrep.preprocess(&mtop_);
 201
 202     // Get data about changes and check it
 203     QMMMTopologyInfo info = topPrep.topInfo();
 204
 205     gmx::test::TestReferenceData    data;
 206     QMMMTopologyPreprocessorChecker checker(data.rootChecker());
 207     checker.checkAll(info, mtop_);
 208 }
 209
 210 TEST_F(QMMMTopologyPreprocessorTest, AlanineDipeptideWithLinksWithConstraints)
 211 {
 212     // Reference input alanine_vacuo.top with constraints=all-bonds
 213     makeMtopFromFile("alanine_vacuo", "constraints = all-bonds");
 214     qmIndices_ = { 8, 9, 10, 11, 12, 13 };
 215
 216     QMMMTopologyPreprocessor topPrep(qmIndices_);
 217     topPrep.preprocess(&mtop_);
 218
 219     // Get data about changes and check it
 220     QMMMTopologyInfo info = topPrep.topInfo();
 221
 222     gmx::test::TestReferenceData    data;
 223     QMMMTopologyPreprocessorChecker checker(data.rootChecker());
 224     checker.checkAll(info, mtop_);
 225 }
 226
 227 TEST_F(QMMMTopologyPreprocessorTest, RemovingQMVsites)
 228 {
 229     // Reference input vistes_test.top
 230     makeMtopFromFile("vsite_test", "");
 231     qmIndices_ = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15 };
 232
 233     QMMMTopologyPreprocessor topPrep(qmIndices_);
 234     topPrep.preprocess(&mtop_);
 235
 236     // Get data about changes and check it
 237     QMMMTopologyInfo info = topPrep.topInfo();
 238
 239     gmx::test::TestReferenceData    data;
 240     QMMMTopologyPreprocessorChecker checker(data.rootChecker());
 241     checker.checkAll(info, mtop_);
 242 }
 243
 244 } // namespace gmx