src/gromacs/gmxana/tests/gmx_msd.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 gmx msd.
  38  *
  39  * \author Kevin Boyd <kevin.boyd@uconn.edu>
  40  */
  41
  42 #include "gmxpre.h"
  43
  44 #include <cstdio>
  45 #include <cstdlib>
  46
  47 #include "gromacs/gmxana/gmx_ana.h"
  48 #include "gromacs/gmxpreprocess/grompp.h"
  49 #include "gromacs/utility/futil.h"
  50 #include "gromacs/utility/path.h"
  51 #include "gromacs/utility/textreader.h"
  52
  53 #include "testutils/cmdlinetest.h"
  54 #include "testutils/refdata.h"
  55 #include "testutils/testfilemanager.h"
  56 #include "testutils/textblockmatchers.h"
  57 #include "testutils/xvgtest.h"
  58
  59 namespace
  60 {
  61
  62 using gmx::test::CommandLine;
  63 using gmx::test::XvgMatch;
  64
  65 class MsdTest : public gmx::test::CommandLineTestBase
  66 {
  67     public:
  68         MsdTest()
  69         {
  70             setOutputFile("-o", "msd.xvg", XvgMatch());
  71             setInputFile("-f", "msd_traj.xtc");
  72             setInputFile("-s", "msd_coords.gro");
  73             setInputFile("-n", "msd.ndx");
  74         }
  75
  76         void runTest(const CommandLine &args)
  77         {
  78             CommandLine &cmdline = commandLine();
  79             cmdline.merge(args);
  80             ASSERT_EQ(0, gmx_msd(cmdline.argc(), cmdline.argv()));
  81             checkOutputFiles();
  82         }
  83 };
  84
  85 class MsdMolTest : public gmx::test::CommandLineTestBase
  86 {
  87     public:
  88         MsdMolTest()
  89         {
  90             double    tolerance = 1e-5;
  91             XvgMatch  xvg;
  92             XvgMatch &toler     = xvg.tolerance(gmx::test::relativeToleranceAsFloatingPoint(1, tolerance));
  93             setOutputFile("-mol", "msdmol.xvg", toler);
  94         }
  95
  96         void runTest(const CommandLine &args, const char *ndxfile,
  97                      const std::string &simulationName)
  98         {
  99             setInputFile("-f", simulationName + ".pdb");
 100             std::string tpr = fileManager().getTemporaryFilePath(".tpr");
 101             std::string mdp = fileManager().getTemporaryFilePath(".mdp");
 102             FILE       *fp  = fopen(mdp.c_str(), "w");
 103             fprintf(fp, "cutoff-scheme = verlet\n");
 104             fprintf(fp, "rcoulomb      = 0.85\n");
 105             fprintf(fp, "rvdw          = 0.85\n");
 106             fprintf(fp, "rlist         = 0.85\n");
 107             fclose(fp);
 108
 109             // Prepare a .tpr file
 110             {
 111                 CommandLine caller;
 112                 auto        simDB = gmx::test::TestFileManager::getTestSimulationDatabaseDirectory();
 113                 auto        base  = gmx::Path::join(simDB, simulationName);
 114                 caller.append("grompp");
 115                 caller.addOption("-maxwarn", 0);
 116                 caller.addOption("-f", mdp.c_str());
 117                 std::string gro = (base + ".pdb");
 118                 caller.addOption("-c", gro.c_str());
 119                 std::string top = (base + ".top");
 120                 caller.addOption("-p", top.c_str());
 121                 std::string ndx = (base + ".ndx");
 122                 caller.addOption("-n", ndx.c_str());
 123                 caller.addOption("-o", tpr.c_str());
 124                 ASSERT_EQ(0, gmx_grompp(caller.argc(), caller.argv()));
 125             }
 126             // Run the MSD analysis
 127             {
 128                 setInputFile("-n", ndxfile);
 129                 CommandLine &cmdline = commandLine();
 130                 cmdline.merge(args);
 131                 cmdline.addOption("-s", tpr.c_str());
 132                 ASSERT_EQ(0, gmx_msd(cmdline.argc(), cmdline.argv()));
 133                 checkOutputFiles();
 134             }
 135         }
 136 };
 137
 138 /* msd_traj.xtc contains a 10 frame (1 ps per frame) simulation
 139  * containing 3 atoms, with different starting positions but identical
 140  * displacements. The displacements are calculated to yield the following
 141  * diffusion coefficients when lag is calculated ONLY FROM TIME 0
 142  * D_x = 8 * 10 ^ -5 cm^2 /s, D_y = 4 * 10^ -5 cm^2 /s , D_z = 0
 143  *
 144  * To test for these results, -trestart is set to a larger value than the
 145  * total simulation length, so that only lag 0 is calculated
 146  */
 147
 148 // for 3D, (8 + 4 + 0) / 3 should yield 4 cm^2 / s
 149 TEST_F(MsdTest, threeDimensionalDiffusion)
 150 {
 151     const char *const cmdline[] = {
 152         "msd", "-mw", "no", "-trestart", "200",
 153     };
 154     runTest(CommandLine(cmdline));
 155 }
 156
 157 // for lateral z, (8 + 4) / 2 should yield 6 cm^2 /s
 158 TEST_F(MsdTest, twoDimensionalDiffusion)
 159 {
 160     const char *const cmdline[] = {
 161         "msd", "-mw", "no", "-trestart", "200", "-lateral", "z"
 162     };
 163     runTest(CommandLine(cmdline));
 164 }
 165
 166 // for type x, should yield 8 cm^2 / s
 167 TEST_F(MsdTest, oneDimensionalDiffusion)
 168 {
 169     const char *const cmdline[] = {
 170         "msd", "-mw", "no", "-trestart", "200", "-type", "x"
 171     };
 172     runTest(CommandLine(cmdline));
 173 }
 174
 175 // Test the diffusion per molecule output, mass weighted
 176 TEST_F(MsdMolTest, diffMolMassWeighted)
 177 {
 178     const char *const cmdline[] = {
 179         "msd", "-trestart", "200"
 180     };
 181     runTest(CommandLine(cmdline), "spc5.ndx", "spc5");
 182 }
 183
 184 // Test the diffusion per molecule output, non-mass weighted
 185 TEST_F(MsdMolTest, diffMolNonMassWeighted)
 186 {
 187     const char *const cmdline[] = {
 188         "msd", "-trestart", "200", "-mw", "no"
 189     };
 190     runTest(CommandLine(cmdline), "spc5.ndx", "spc5");
 191 }
 192
 193 // Test the diffusion per molecule output, with selection
 194 TEST_F(MsdMolTest, diffMolSelected)
 195 {
 196     const char *const cmdline[] = {
 197         "msd", "-trestart", "200"
 198     };
 199     runTest(CommandLine(cmdline), "spc5_3.ndx", "spc5");
 200 }
 201
 202 } //namespace