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, const std::string& simulationName)
  97     {
  98         setInputFile("-f", simulationName + ".pdb");
  99         std::string tpr = fileManager().getTemporaryFilePath(".tpr");
 100         std::string mdp = fileManager().getTemporaryFilePath(".mdp");
 101         FILE*       fp  = fopen(mdp.c_str(), "w");
 102         fprintf(fp, "cutoff-scheme = verlet\n");
 103         fprintf(fp, "rcoulomb      = 0.85\n");
 104         fprintf(fp, "rvdw          = 0.85\n");
 105         fprintf(fp, "rlist         = 0.85\n");
 106         fclose(fp);
 107
 108         // Prepare a .tpr file
 109         {
 110             CommandLine caller;
 111             auto        simDB = gmx::test::TestFileManager::getTestSimulationDatabaseDirectory();
 112             auto        base  = gmx::Path::join(simDB, simulationName);
 113             caller.append("grompp");
 114             caller.addOption("-maxwarn", 0);
 115             caller.addOption("-f", mdp.c_str());
 116             std::string gro = (base + ".pdb");
 117             caller.addOption("-c", gro.c_str());
 118             std::string top = (base + ".top");
 119             caller.addOption("-p", top.c_str());
 120             std::string ndx = (base + ".ndx");
 121             caller.addOption("-n", ndx.c_str());
 122             caller.addOption("-o", tpr.c_str());
 123             ASSERT_EQ(0, gmx_grompp(caller.argc(), caller.argv()));
 124         }
 125         // Run the MSD analysis
 126         {
 127             setInputFile("-n", ndxfile);
 128             CommandLine& cmdline = commandLine();
 129             cmdline.merge(args);
 130             cmdline.addOption("-s", tpr.c_str());
 131             ASSERT_EQ(0, gmx_msd(cmdline.argc(), cmdline.argv()));
 132             checkOutputFiles();
 133         }
 134     }
 135 };
 136
 137 /* msd_traj.xtc contains a 10 frame (1 ps per frame) simulation
 138  * containing 3 atoms, with different starting positions but identical
 139  * displacements. The displacements are calculated to yield the following
 140  * diffusion coefficients when lag is calculated ONLY FROM TIME 0
 141  * D_x = 8 * 10 ^ -5 cm^2 /s, D_y = 4 * 10^ -5 cm^2 /s , D_z = 0
 142  *
 143  * To test for these results, -trestart is set to a larger value than the
 144  * total simulation length, so that only lag 0 is calculated
 145  */
 146
 147 // for 3D, (8 + 4 + 0) / 3 should yield 4 cm^2 / s
 148 TEST_F(MsdTest, threeDimensionalDiffusion)
 149 {
 150     const char* const cmdline[] = {
 151         "msd", "-mw", "no", "-trestart", "200",
 152     };
 153     runTest(CommandLine(cmdline));
 154 }
 155
 156 // for lateral z, (8 + 4) / 2 should yield 6 cm^2 /s
 157 TEST_F(MsdTest, twoDimensionalDiffusion)
 158 {
 159     const char* const cmdline[] = { "msd", "-mw", "no", "-trestart", "200", "-lateral", "z" };
 160     runTest(CommandLine(cmdline));
 161 }
 162
 163 // for type x, should yield 8 cm^2 / s
 164 TEST_F(MsdTest, oneDimensionalDiffusion)
 165 {
 166     const char* const cmdline[] = { "msd", "-mw", "no", "-trestart", "200", "-type", "x" };
 167     runTest(CommandLine(cmdline));
 168 }
 169
 170 // Test the diffusion per molecule output, mass weighted
 171 TEST_F(MsdMolTest, diffMolMassWeighted)
 172 {
 173     const char* const cmdline[] = { "msd", "-trestart", "200" };
 174     runTest(CommandLine(cmdline), "spc5.ndx", "spc5");
 175 }
 176
 177 // Test the diffusion per molecule output, non-mass weighted
 178 TEST_F(MsdMolTest, diffMolNonMassWeighted)
 179 {
 180     const char* const cmdline[] = { "msd", "-trestart", "200", "-mw", "no" };
 181     runTest(CommandLine(cmdline), "spc5.ndx", "spc5");
 182 }
 183
 184 // Test the diffusion per molecule output, with selection
 185 TEST_F(MsdMolTest, diffMolSelected)
 186 {
 187     const char* const cmdline[] = { "msd", "-trestart", "200" };
 188     runTest(CommandLine(cmdline), "spc5_3.ndx", "spc5");
 189 }
 190
 191 } // namespace