Add PME integration test
[alexxy/gromacs.git] / src / programs / mdrun / tests / pmetest.cpp
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
4  * Copyright (c) 2016,2017, 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  * This implements basic PME sanity test (using single-rank mdrun).
38  * It runs the input system with PME for several steps (on CPU and GPU, if available),
39  * and checks the reciprocal and conserved energies.
40  * TODO: implement multi-rank tests as well.
41  *
42  * \author Aleksei Iupinov <a.yupinov@gmail.com>
43  * \ingroup module_mdrun_integration_tests
44  */
45 #include "gmxpre.h"
46
47 #include "config.h"
48
49 #include <map>
50 #include <string>
51 #include <vector>
52
53 #include "gromacs/gpu_utils/gpu_utils.h"
54 #include "gromacs/hardware/gpu_hw_info.h"
55 #include "gromacs/utility/cstringutil.h"
56 #include "gromacs/utility/stringutil.h"
57
58 #include "testutils/refdata.h"
59
60 #include "energyreader.h"
61 #include "moduletest.h"
62
63 namespace gmx
64 {
65 namespace test
66 {
67 namespace
68 {
69
70 //! A basic PME runner
71 class PmeTest : public MdrunTestFixture
72 {
73     public:
74         //! Before any test is run, work out whether any compatible GPUs exist.
75         static void SetUpTestCase();
76         //! Store whether any compatible GPUs exist.
77         static bool s_hasCompatibleCudaGpus;
78 };
79
80 bool PmeTest::s_hasCompatibleCudaGpus = false;
81
82 void PmeTest::SetUpTestCase()
83 {
84     gmx_gpu_info_t gpuInfo {};
85     char           detection_error[STRLEN];
86     GMX_UNUSED_VALUE(detection_error); //TODO
87     // It would be nicer to do this detection once and have mdrun
88     // re-use it, but this is OK. Note that this also caters for when
89     // there is no GPU support in the build.
90     if (GMX_GPU == GMX_GPU_CUDA &&
91         (detect_gpus(&gpuInfo, detection_error) >= 0) &&
92         gpuInfo.n_dev_compatible > 0)
93     {
94         s_hasCompatibleCudaGpus = true;
95     }
96     free_gpu_info(&gpuInfo);
97 }
98
99 TEST_F(PmeTest, ReproducesEnergies)
100 {
101     const int   nsteps     = 20;
102     std::string theMdpFile = formatString("coulombtype     = PME\n"
103                                           "nstcalcenergy   = 1\n"
104                                           "nstenergy       = 1\n"
105                                           "pme-order       = 4\n"
106                                           "nsteps          = %d\n",
107                                           nsteps);
108
109     runner_.useStringAsMdpFile(theMdpFile);
110
111     const std::string inputFile = "spc-and-methanol";
112     runner_.useTopGroAndNdxFromDatabase(inputFile.c_str());
113
114     EXPECT_EQ(0, runner_.callGrompp());
115
116     //TODO test all proper/improper combinations in more thorough way?
117     std::map < std::string, std::vector < const char *>> runModes;
118     runModes["PmeOnCpu"] = {};
119     // TODO uncomment these and replace the above as functionality
120     // gets implemented.
121     //runModes["PmeOnCpu"] = {"-pme", "cpu"};
122     //runModes["PmeAuto"]  = {"-pme", "auto"};
123     //runModes["PmeOnGpuFftOnCpu"] = {"-pme", "gpu", "-pmefft", "cpu"};
124     //runModes["PmeOnGpuFftOnGpu"] = {"-pme", "gpu", "-pmefft", "gpu"};
125     //runModes["PmeOnGpuFftAuto"] = {"-pme", "gpu", "-pmefft", "auto"};
126     TestReferenceData    refData;
127     TestReferenceChecker rootChecker(refData.rootChecker());
128
129     for (const auto &mode : runModes)
130     {
131         auto modeTargetsGpus = (mode.first.find("Gpu") != std::string::npos);
132         if (modeTargetsGpus && !s_hasCompatibleCudaGpus)
133         {
134             // This run mode will cause a fatal error from mdrun when
135             // it can't find GPUs, which is not something we're trying
136             // to test here.
137             continue;
138         }
139
140         runner_.edrFileName_ = fileManager_.getTemporaryFilePath(inputFile + "_" + mode.first + ".edr");
141
142         CommandLine commandLine(mode.second);
143         commandLine.append("-notunepme"); // for reciprocal energy reproducibility
144         ASSERT_EQ(0, runner_.callMdrun(commandLine));
145
146         auto energyReader      = openEnergyFileToReadFields(runner_.edrFileName_, {"Coul. recip.", "Total Energy", "Kinetic En."});
147         auto conservedChecker  = rootChecker.checkCompound("Energy", "Conserved");
148         auto reciprocalChecker = rootChecker.checkCompound("Energy", "Reciprocal");
149         for (int i = 0; i <= nsteps; i++)
150         {
151             EnergyFrame frame            = energyReader->frame();
152             std::string stepNum          = gmx::formatString("%d", i);
153             const real  conservedEnergy  = frame.at("Total Energy");
154             const real  reciprocalEnergy = frame.at("Coul. recip.");
155             if (i == 0)
156             {
157                 // use first step values as references for tolerance
158                 const real startingKineticEnergy = frame.at("Kinetic En.");
159                 const auto conservedTolerance    = relativeToleranceAsFloatingPoint(startingKineticEnergy, 2e-5);
160                 const auto reciprocalTolerance   = relativeToleranceAsFloatingPoint(reciprocalEnergy, 2e-5);
161                 reciprocalChecker.setDefaultTolerance(reciprocalTolerance);
162                 conservedChecker.setDefaultTolerance(conservedTolerance);
163             }
164             conservedChecker.checkReal(conservedEnergy, stepNum.c_str());
165             reciprocalChecker.checkReal(reciprocalEnergy, stepNum.c_str());
166         }
167     }
168 }
169
170 }
171 }
172 }