src/gromacs/mdlib/tests/energyoutput.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 energy output to log and .edr files.
  38  *
  39  * \todo Position and orientation restraints tests.
  40  * \todo Average and sum in edr file test.
  41  * \todo AWH output tests.
  42  * \todo The log and edr outputs are currently saved to the file on the disk and then read
  43  *       to compare with the reference data. This will be more elegant (and run faster) when we
  44  *       refactor the output routines to write to a stream interface, which can already be handled
  45  *       in-memory when running tests.
  46  *
  47  * \author Mark Abraham <mark.j.abraham@gmail.com>
  48  * \author Artem Zhmurov <zhmurov@gmail.com>
  49  *
  50  * \ingroup module_mdlib
  51  */
  52 #include "gmxpre.h"
  53
  54 #include "gromacs/mdlib/energyoutput.h"
  55
  56 #include "config.h"
  57
  58 #include <cstdio>
  59 #include <cstdlib>
  60
  61 #include <gtest/gtest.h>
  62
  63 #include "gromacs/mdlib/ebin.h"
  64 #include "gromacs/mdlib/makeconstraints.h"
  65 #include "gromacs/mdtypes/commrec.h"
  66 #include "gromacs/mdtypes/fcdata.h"
  67 #include "gromacs/mdtypes/group.h"
  68 #include "gromacs/mdtypes/inputrec.h"
  69 #include "gromacs/mdtypes/mdatom.h"
  70 #include "gromacs/mdtypes/state.h"
  71 #include "gromacs/topology/topology.h"
  72 #include "gromacs/utility/cstringutil.h"
  73 #include "gromacs/utility/stringutil.h"
  74 #include "gromacs/utility/textreader.h"
  75 #include "gromacs/utility/unique_cptr.h"
  76
  77 #include "testutils/refdata.h"
  78 #include "testutils/testasserts.h"
  79 #include "testutils/testfilemanager.h"
  80
  81 namespace gmx
  82 {
  83 namespace test
  84 {
  85 namespace
  86 {
  87
  88 //! Wraps fclose to discard the return value to use it as a deleter with gmx::unique_cptr.
  89 void fcloseWrapper(FILE *fp)
  90 {
  91     fclose(fp);
  92 }
  93
  94 /*! \brief Test parameters space.
  95  *
  96  * The test will run on a set of combinations of this steucture parameters.
  97  */
  98 struct EnergyOutputTestParameters
  99 {
 100     //! If output should be initialized as a rerun.
 101     bool isRerun;
 102     //! Thermostat (enum)
 103     int  temperatureCouplingScheme;
 104     //! Barostat (enum)
 105     int  pressureCouplingScheme;
 106     //! Integrator
 107     int  integrator;
 108     //! Is box triclinic (off-diagonal elements will be printed).
 109     bool isBoxTriclinic;
 110     //! Number of saved energy frames (to test averages output).
 111     int  numFrames;
 112 };
 113
 114 /*! \brief Sets of parameters on which to run the tests.
 115  *
 116  * Only several combinations of the parameters are used. Using all possible combinations will require ~10 MB of
 117  * test data and ~2 sec to run the tests.
 118  */
 119 const EnergyOutputTestParameters parametersSets[] = {{false, etcNO,         epcNO,               eiMD, false, 1},
 120                                                      {true,  etcNO,         epcNO,               eiMD, false, 1},
 121                                                      {false, etcNO,         epcNO,               eiMD, true,  1},
 122                                                      {false, etcNO,         epcNO,               eiMD, false, 0},
 123                                                      {false, etcNO,         epcNO,               eiMD, false, 10},
 124                                                      {false, etcVRESCALE,   epcNO,               eiMD, false, 1},
 125                                                      {false, etcNOSEHOOVER, epcNO,               eiMD, false, 1},
 126                                                      {false, etcNO,         epcPARRINELLORAHMAN, eiMD, false, 1},
 127                                                      {false, etcNO,         epcMTTK,             eiMD, false, 1},
 128                                                      {false, etcNO,         epcNO,               eiVV, false, 1},
 129                                                      {false, etcNO,         epcMTTK,             eiVV, false, 1}};
 130
 131 /*! \brief Test fixture to test energy output.
 132  *
 133  * The class is initialized to maximize amount of energy terms printed.
 134  * The test is run for different combinations of temperature and pressure control
 135  * schemes. Different number of printed steps is also tested.
 136  */
 137 class EnergyOutputTest : public ::testing::TestWithParam<EnergyOutputTestParameters>
 138 {
 139     public:
 140         //! File manager
 141         TestFileManager                     fileManager_;
 142         //! Energy (.edr) file
 143         ener_file_t                         energyFile_;
 144         //! Input data
 145         t_inputrec                          inputrec_;
 146         //! Topology
 147         gmx_mtop_t                          mtop_;
 148         //! Energy output object
 149         EnergyOutput                        energyOutput_;
 150         //! MD atoms
 151         t_mdatoms                           mdatoms_;
 152         //! Simulation time
 153         double                              time_;
 154         //! Total mass
 155         real                                tmass_;
 156         //! Potential energy data
 157         std::unique_ptr<gmx_enerdata_t>     enerdata_;
 158         //! Kinetic energy data (for temperatures output)
 159         gmx_ekindata_t                      ekindata_;
 160         //! System state
 161         t_state                             state_;
 162         //! PBC box
 163         matrix                              box_;
 164         //! Virial from constraints
 165         tensor                              constraintsVirial_;
 166         //! Virial from force computation
 167         tensor                              forceVirial_;
 168         //! Total virial
 169         tensor                              totalVirial_;
 170         //! Pressure
 171         tensor                              pressure_;
 172         //! Names for the groups.
 173         std::vector<std::string>            groupNameStrings_ = { "Protein", "Water", "Lipid" };
 174         //! Names for the groups as C strings.
 175         std::vector < std::vector < char>>  groupNameCStrings_;
 176         //! Handles to the names as C strings in the way needed for SimulationGroups.
 177         std::vector<char *>                 groupNameHandles_;
 178         //! Total dipole momentum
 179         rvec                                muTotal_;
 180         //! Distance and orientation restraints data
 181         t_fcdata                            fcd_;
 182         //! Communication record
 183         t_commrec                           cr_;
 184         //! Constraints object (for constraints RMSD output in case of LINCS)
 185         std::unique_ptr<Constraints>        constraints_;
 186         //! \brief Temporary output filename
 187         std::string                         logFilename_;
 188         //! Temporary energy output filename
 189         std::string                         edrFilename_;
 190         //! Pointer to a temporary output file
 191         FILE                               *log_;
 192         //! Log file wrapper
 193         unique_cptr<FILE, fcloseWrapper>    logFileGuard_;
 194         //! Reference data
 195         TestReferenceData                   refData_;
 196         //! Checker for reference data
 197         TestReferenceChecker                checker_;
 198
 199         EnergyOutputTest() :
 200             logFilename_(fileManager_.getTemporaryFilePath(".log")),
 201             edrFilename_(fileManager_.getTemporaryFilePath(".edr")),
 202             log_(std::fopen(logFilename_.c_str(), "w")),
 203             logFileGuard_(log_),
 204             checker_(refData_.rootChecker())
 205         {
 206             // Input record
 207             inputrec_.delta_t = 0.001;
 208
 209             // F_EQM
 210             inputrec_.bQMMM          = true;
 211             // F_RF_EXCL will not be tested - group scheme is not supported any more
 212             inputrec_.cutoff_scheme  = ecutsVERLET;
 213             // F_COUL_RECIP
 214             inputrec_.coulombtype    = eelPME;
 215             // F_LJ_RECIP
 216             inputrec_.vdwtype        = evdwPME;
 217
 218             // F_DVDL_COUL, F_DVDL_VDW, F_DVDL_BONDED, F_DVDL_RESTRAINT, F_DKDL and F_DVDL
 219             inputrec_.efep = efepYES;
 220             inputrec_.fepvals->separate_dvdl[efptCOUL]      = true;
 221             inputrec_.fepvals->separate_dvdl[efptVDW]       = true;
 222             inputrec_.fepvals->separate_dvdl[efptBONDED]    = true;
 223             inputrec_.fepvals->separate_dvdl[efptRESTRAINT] = true;
 224             inputrec_.fepvals->separate_dvdl[efptMASS]      = true;
 225             inputrec_.fepvals->separate_dvdl[efptCOUL]      = true;
 226             inputrec_.fepvals->separate_dvdl[efptFEP]       = true;
 227
 228             // F_DISPCORR and F_PDISPCORR
 229             inputrec_.eDispCorr       = edispcEner;
 230             inputrec_.bRot            = true;
 231
 232             // F_ECONSERVED
 233             inputrec_.ref_p[YY][XX]   = 0.0;
 234             inputrec_.ref_p[ZZ][XX]   = 0.0;
 235             inputrec_.ref_p[ZZ][YY]   = 0.0;
 236
 237             // Dipole (mu)
 238             inputrec_.ewald_geometry  = eewg3DC;
 239
 240             // GMX_CONSTRAINTVIR environment variable should also be
 241             // set to print constraints and force virials separately.
 242             //
 243             // TODO extract a helper function if we ever need to do
 244             // this kind of thing again.
 245 #if GMX_NATIVE_WINDOWS
 246             _putenv_s("GMX_CONSTRAINTVIR", "true");
 247 #else
 248             setenv("GMX_CONSTRAINTVIR", "true", 1);
 249 #endif
 250             // To print constrain RMSD, constraints algorithm should be set to LINCS.
 251             inputrec_.eConstrAlg = econtLINCS;
 252
 253             mtop_.bIntermolecularInteractions = false;
 254
 255             // Constructing molecular topology
 256             gmx_moltype_t molType;
 257
 258             molType.atoms.nr = 2;
 259
 260             // F_CONSTR
 261             // This must be initialized so that Constraints object can be created below.
 262             InteractionList interactionListConstr;
 263             interactionListConstr.iatoms.resize(NRAL(F_CONSTR) + 1);
 264             interactionListConstr.iatoms[0] = 0;
 265             interactionListConstr.iatoms[1] = 0;
 266             interactionListConstr.iatoms[2] = 1;
 267             molType.ilist.at(F_CONSTR)      = interactionListConstr;
 268
 269             InteractionList interactionListEmpty;
 270             interactionListEmpty.iatoms.resize(0);
 271             molType.ilist.at(F_CONSTRNC)   = interactionListEmpty;
 272             molType.ilist.at(F_SETTLE)     = interactionListEmpty;
 273
 274             // F_LJ14 and F_COUL14
 275             InteractionList interactionListLJ14;
 276             interactionListLJ14.iatoms.resize(NRAL(F_LJ14) + 1);
 277             molType.ilist.at(F_LJ14)     = interactionListLJ14;
 278
 279             // F_LJC14_Q
 280             InteractionList interactionListLJC14Q;
 281             interactionListLJC14Q.iatoms.resize(NRAL(F_LJC14_Q) + 1);
 282             molType.ilist.at(F_LJC14_Q)  = interactionListLJC14Q;
 283
 284             // TODO Do proper initialization for distance and orientation
 285             //      restraints and remove comments to enable their output
 286             // F_DISRES
 287             //InteractionList interactionListDISRES;
 288             //interactionListDISRES.iatoms.resize(NRAL(F_DISRES) + 1);
 289             //molType.ilist.at(F_DISRES)   = interactionListDISRES;
 290             //
 291             // F_ORIRES
 292             //InteractionList interactionListORIRES;
 293             //interactionListORIRES.iatoms.resize(NRAL(F_ORIRES) + 1);
 294             //molType.ilist.at(F_ORIRES)   = interactionListORIRES;
 295
 296             mtop_.moltype.push_back(molType);
 297
 298             gmx_molblock_t molBlock;
 299             molBlock.type = 0;
 300             molBlock.nmol = 1;
 301             mtop_.molblock.push_back(molBlock);
 302
 303             // This is to keep constraints initialization happy
 304             mtop_.natoms = 2;
 305             mtop_.ffparams.iparams.resize(F_NRE);
 306             mtop_.ffparams.functype.resize(F_NRE);
 307             mtop_.ffparams.iparams.at(F_CONSTR).constr.dA   = 1.0;
 308             mtop_.ffparams.iparams.at(F_CONSTR).constr.dB   = 1.0;
 309             mtop_.ffparams.iparams.at(F_CONSTRNC).constr.dA = 1.0;
 310             mtop_.ffparams.iparams.at(F_CONSTRNC).constr.dB = 1.0;
 311
 312             // Groups for energy output, temperature coupling and acceleration
 313             for (const auto &string : groupNameStrings_)
 314             {
 315                 std::vector<char> cString(string.begin(), string.end());
 316                 // Need to add null termination
 317                 cString.push_back('\0');
 318                 groupNameCStrings_.emplace_back(cString);
 319                 groupNameHandles_.emplace_back(groupNameCStrings_.back().data());
 320             }
 321             for (auto &handle : groupNameHandles_)
 322             {
 323                 mtop_.groups.groupNames.emplace_back(&handle);
 324             }
 325
 326             mtop_.groups.groups[SimulationAtomGroupType::EnergyOutput].nr = 3;
 327             snew(mtop_.groups.groups[SimulationAtomGroupType::EnergyOutput].nm_ind,
 328                  mtop_.groups.groups[SimulationAtomGroupType::EnergyOutput].nr);
 329             mtop_.groups.groups[SimulationAtomGroupType::EnergyOutput].nm_ind[0] = 0;
 330             mtop_.groups.groups[SimulationAtomGroupType::EnergyOutput].nm_ind[1] = 1;
 331             mtop_.groups.groups[SimulationAtomGroupType::EnergyOutput].nm_ind[2] = 2;
 332
 333             mtop_.groups.groups[SimulationAtomGroupType::TemperatureCoupling].nr = 3;
 334             snew(mtop_.groups.groups[SimulationAtomGroupType::TemperatureCoupling].nm_ind,
 335                  mtop_.groups.groups[SimulationAtomGroupType::TemperatureCoupling].nr);
 336             mtop_.groups.groups[SimulationAtomGroupType::TemperatureCoupling].nm_ind[0] = 0;
 337             mtop_.groups.groups[SimulationAtomGroupType::TemperatureCoupling].nm_ind[1] = 1;
 338             mtop_.groups.groups[SimulationAtomGroupType::TemperatureCoupling].nm_ind[2] = 2;
 339
 340             mtop_.groups.groups[SimulationAtomGroupType::Acceleration].nr = 2;
 341             snew(mtop_.groups.groups[SimulationAtomGroupType::Acceleration].nm_ind,
 342                  mtop_.groups.groups[SimulationAtomGroupType::Acceleration].nr);
 343             mtop_.groups.groups[SimulationAtomGroupType::Acceleration].nm_ind[0] = 0;
 344             mtop_.groups.groups[SimulationAtomGroupType::Acceleration].nm_ind[1] = 2;
 345
 346             // Nose-Hoover chains
 347             inputrec_.bPrintNHChains     = true;
 348             inputrec_.opts.nhchainlength = 2;
 349             state_.nosehoover_xi.resize(mtop_.groups.groups[SimulationAtomGroupType::TemperatureCoupling].nr*inputrec_.opts.nhchainlength);
 350             state_.nosehoover_vxi.resize(mtop_.groups.groups[SimulationAtomGroupType::TemperatureCoupling].nr*inputrec_.opts.nhchainlength);
 351
 352             // This will be needed only with MTTK barostat
 353             state_.nhpres_xi.resize(1*inputrec_.opts.nhchainlength);
 354             state_.nhpres_vxi.resize(1*inputrec_.opts.nhchainlength);
 355
 356             // Group pairs
 357             enerdata_ = std::make_unique<gmx_enerdata_t>(mtop_.groups.groups[SimulationAtomGroupType::EnergyOutput].nr, 0);
 358
 359             // Kinetic energy and related data
 360             ekindata_.tcstat.resize(mtop_.groups.groups[SimulationAtomGroupType::TemperatureCoupling].nr);
 361             ekindata_.grpstat.resize(mtop_.groups.groups[SimulationAtomGroupType::Acceleration].nr);
 362
 363             // This is needed so that the ebin space will be allocated
 364             inputrec_.cos_accel = 1.0;
 365             // This is to keep the destructor happy (otherwise sfree() segfaults)
 366             ekindata_.nthreads = 0;
 367             snew(ekindata_.ekin_work_alloc, 1);
 368             snew(ekindata_.ekin_work, 1);
 369             snew(ekindata_.dekindl_work, 1);
 370
 371             // Group options for annealing output
 372             inputrec_.opts.ngtc = 3;
 373             snew(inputrec_.opts.ref_t, inputrec_.opts.ngtc);
 374             snew(inputrec_.opts.annealing, inputrec_.opts.ngtc);
 375             inputrec_.opts.annealing[0] = eannNO;
 376             inputrec_.opts.annealing[1] = eannSINGLE;
 377             inputrec_.opts.annealing[2] = eannPERIODIC;
 378
 379             // This is to keep done_inputrec happy (otherwise sfree() segfaults)
 380             snew(inputrec_.opts.anneal_time, inputrec_.opts.ngtc);
 381             snew(inputrec_.opts.anneal_temp, inputrec_.opts.ngtc);
 382
 383             // Communication record (for Constraints constructor)
 384             cr_.nnodes = 1;
 385             cr_.dd     = nullptr;
 386
 387             // Constraints object (to get constraints RMSD from)
 388             // TODO EnergyOutput should not take Constraints object
 389             // TODO This object will always return zero as RMSD value.
 390             //      It is more relevant to have non-zero value for testing.
 391             constraints_ = makeConstraints(mtop_, inputrec_, nullptr, false, nullptr, mdatoms_, &cr_,
 392                                            nullptr, nullptr, nullptr, false);
 393
 394             // No position/orientation restraints
 395             fcd_.disres.npair  = 0;
 396             fcd_.orires.nr     = 0;
 397
 398         }
 399
 400         /*! \brief Helper function to generate synthetic data to output
 401          *
 402          * \param[in,out] testValue    Base value fr energy data.
 403          */
 404         void setStepData(double *testValue)
 405         {
 406
 407             time_                             = (*testValue += 0.1);
 408             tmass_                            = (*testValue += 0.1);
 409
 410             enerdata_->term[F_LJ]             = (*testValue += 0.1);
 411             enerdata_->term[F_COUL_SR]        = (*testValue += 0.1);
 412             enerdata_->term[F_EPOT]           = (*testValue += 0.1);
 413             enerdata_->term[F_EKIN]           = (*testValue += 0.1);
 414             enerdata_->term[F_ETOT]           = (*testValue += 0.1);
 415             enerdata_->term[F_TEMP]           = (*testValue += 0.1);
 416             enerdata_->term[F_PRES]           = (*testValue += 0.1);
 417
 418             enerdata_->term[F_BHAM]           = (*testValue += 0.1);
 419             enerdata_->term[F_EQM]            = (*testValue += 0.1);
 420             enerdata_->term[F_RF_EXCL]        = (*testValue += 0.1);
 421             enerdata_->term[F_COUL_RECIP]     = (*testValue += 0.1);
 422             enerdata_->term[F_LJ_RECIP]       = (*testValue += 0.1);
 423             enerdata_->term[F_LJ14]           = (*testValue += 0.1);
 424             enerdata_->term[F_COUL14]         = (*testValue += 0.1);
 425             enerdata_->term[F_LJC14_Q]        = (*testValue += 0.1);
 426             enerdata_->term[F_LJC_PAIRS_NB]   = (*testValue += 0.1);
 427
 428             enerdata_->term[F_DVDL_COUL]      = (*testValue += 0.1);
 429             enerdata_->term[F_DVDL_VDW]       = (*testValue += 0.1);
 430             enerdata_->term[F_DVDL_BONDED]    = (*testValue += 0.1);
 431             enerdata_->term[F_DVDL_RESTRAINT] = (*testValue += 0.1);
 432             enerdata_->term[F_DKDL]           = (*testValue += 0.1);
 433             enerdata_->term[F_DVDL]           = (*testValue += 0.1);
 434
 435             enerdata_->term[F_DISPCORR]       = (*testValue += 0.1);
 436             enerdata_->term[F_PDISPCORR]      = (*testValue += 0.1);
 437             enerdata_->term[F_DISRESVIOL]     = (*testValue += 0.1);
 438             enerdata_->term[F_ORIRESDEV]      = (*testValue += 0.1);
 439             enerdata_->term[F_COM_PULL]       = (*testValue += 0.1);
 440             enerdata_->term[F_ECONSERVED]     = (*testValue += 0.1);
 441
 442             // Group pairs
 443             for (int i = 0; i < enerdata_->grpp.nener; i++)
 444             {
 445                 for (int k = 0; k < egNR; k++)
 446                 {
 447                     enerdata_->grpp.ener[k][i] = (*testValue += 0.1);
 448                 }
 449             }
 450
 451             // Kinetic energy and related data
 452             for (int i = 0; i < mtop_.groups.groups[SimulationAtomGroupType::TemperatureCoupling].nr; i++)
 453             {
 454                 ekindata_.tcstat[i].T      = (*testValue += 0.1);
 455                 ekindata_.tcstat[i].lambda = (*testValue += 0.1);
 456             }
 457             for (int i = 0; i < mtop_.groups.groups[SimulationAtomGroupType::Acceleration].nr; i++)
 458             {
 459                 ekindata_.grpstat[i].u[XX] = (*testValue += 0.1);
 460                 ekindata_.grpstat[i].u[YY] = (*testValue += 0.1);
 461                 ekindata_.grpstat[i].u[ZZ] = (*testValue += 0.1);
 462             }
 463
 464             // This conditional is to check whether the ebin was allocated.
 465             // Otherwise it will print cosacc data into the first bin.
 466             if (inputrec_.cos_accel != 0)
 467             {
 468                 ekindata_.cosacc.cos_accel = (*testValue += 0.1);
 469                 ekindata_.cosacc.vcos      = (*testValue += 0.1);
 470             }
 471
 472             state_.box[XX][XX]         = (*testValue += 0.1);
 473             state_.box[XX][YY]         = (*testValue += 0.1);
 474             state_.box[XX][ZZ]         = (*testValue += 0.1);
 475             state_.box[YY][XX]         = (*testValue += 0.1);
 476             state_.box[YY][YY]         = (*testValue += 0.1);
 477             state_.box[YY][ZZ]         = (*testValue += 0.1);
 478             state_.box[ZZ][XX]         = (*testValue += 0.1);
 479             state_.box[ZZ][YY]         = (*testValue += 0.1);
 480             state_.box[ZZ][ZZ]         = (*testValue += 0.1);
 481
 482             box_[XX][XX]               = (*testValue += 0.1);
 483             box_[XX][YY]               = (*testValue += 0.1);
 484             box_[XX][ZZ]               = (*testValue += 0.1);
 485             box_[YY][XX]               = (*testValue += 0.1);
 486             box_[YY][YY]               = (*testValue += 0.1);
 487             box_[YY][ZZ]               = (*testValue += 0.1);
 488             box_[ZZ][XX]               = (*testValue += 0.1);
 489             box_[ZZ][YY]               = (*testValue += 0.1);
 490             box_[ZZ][ZZ]               = (*testValue += 0.1);
 491
 492             constraintsVirial_[XX][XX] = (*testValue += 0.1);
 493             constraintsVirial_[XX][YY] = (*testValue += 0.1);
 494             constraintsVirial_[XX][ZZ] = (*testValue += 0.1);
 495             constraintsVirial_[YY][XX] = (*testValue += 0.1);
 496             constraintsVirial_[YY][YY] = (*testValue += 0.1);
 497             constraintsVirial_[YY][ZZ] = (*testValue += 0.1);
 498             constraintsVirial_[ZZ][XX] = (*testValue += 0.1);
 499             constraintsVirial_[ZZ][YY] = (*testValue += 0.1);
 500             constraintsVirial_[ZZ][ZZ] = (*testValue += 0.1);
 501
 502             forceVirial_[XX][XX]       = (*testValue += 0.1);
 503             forceVirial_[XX][YY]       = (*testValue += 0.1);
 504             forceVirial_[XX][ZZ]       = (*testValue += 0.1);
 505             forceVirial_[YY][XX]       = (*testValue += 0.1);
 506             forceVirial_[YY][YY]       = (*testValue += 0.1);
 507             forceVirial_[YY][ZZ]       = (*testValue += 0.1);
 508             forceVirial_[ZZ][XX]       = (*testValue += 0.1);
 509             forceVirial_[ZZ][YY]       = (*testValue += 0.1);
 510             forceVirial_[ZZ][ZZ]       = (*testValue += 0.1);
 511
 512             totalVirial_[XX][XX]       = (*testValue += 0.1);
 513             totalVirial_[XX][YY]       = (*testValue += 0.1);
 514             totalVirial_[XX][ZZ]       = (*testValue += 0.1);
 515             totalVirial_[YY][XX]       = (*testValue += 0.1);
 516             totalVirial_[YY][YY]       = (*testValue += 0.1);
 517             totalVirial_[YY][ZZ]       = (*testValue += 0.1);
 518             totalVirial_[ZZ][XX]       = (*testValue += 0.1);
 519             totalVirial_[ZZ][YY]       = (*testValue += 0.1);
 520             totalVirial_[ZZ][ZZ]       = (*testValue += 0.1);
 521
 522             pressure_[XX][XX]          = (*testValue += 0.1);
 523             pressure_[XX][YY]          = (*testValue += 0.1);
 524             pressure_[XX][ZZ]          = (*testValue += 0.1);
 525             pressure_[YY][XX]          = (*testValue += 0.1);
 526             pressure_[YY][YY]          = (*testValue += 0.1);
 527             pressure_[YY][ZZ]          = (*testValue += 0.1);
 528             pressure_[ZZ][XX]          = (*testValue += 0.1);
 529             pressure_[ZZ][YY]          = (*testValue += 0.1);
 530             pressure_[ZZ][ZZ]          = (*testValue += 0.1);
 531
 532             muTotal_[XX]               = (*testValue += 0.1);
 533             muTotal_[YY]               = (*testValue += 0.1);
 534             muTotal_[ZZ]               = (*testValue += 0.1);
 535
 536             state_.boxv[XX][XX]        = (*testValue += 0.1);
 537             state_.boxv[XX][YY]        = (*testValue += 0.1);
 538             state_.boxv[XX][ZZ]        = (*testValue += 0.1);
 539             state_.boxv[YY][XX]        = (*testValue += 0.1);
 540             state_.boxv[YY][YY]        = (*testValue += 0.1);
 541             state_.boxv[YY][ZZ]        = (*testValue += 0.1);
 542             state_.boxv[ZZ][XX]        = (*testValue += 0.1);
 543             state_.boxv[ZZ][YY]        = (*testValue += 0.1);
 544             state_.boxv[ZZ][ZZ]        = (*testValue += 0.1);
 545
 546             for (int i = 0; i < inputrec_.opts.ngtc; i++)
 547             {
 548                 inputrec_.opts.ref_t[i] = (*testValue += 0.1);
 549             }
 550
 551             for (int k = 0; k < mtop_.groups.groups[SimulationAtomGroupType::TemperatureCoupling].nr*inputrec_.opts.nhchainlength; k++)
 552             {
 553                 state_.nosehoover_xi[k]  = (*testValue += 0.1);
 554                 state_.nosehoover_vxi[k] = (*testValue += 0.1);
 555             }
 556             for (int k = 0; k < inputrec_.opts.nhchainlength; k++)
 557             {
 558                 state_.nhpres_xi[k]  = (*testValue += 0.1);
 559                 state_.nhpres_vxi[k] = (*testValue += 0.1);
 560             }
 561         }
 562
 563         /*! \brief Check if the contents of the .edr file correspond to the reference data.
 564          *
 565          * The code below is based on the 'gmx dump' tool.
 566          *
 567          * \param[in] fileName    Name of the file to check.
 568          * \param[in] frameCount  Number of frames to check.
 569          */
 570         void checkEdrFile(const char *fileName, int frameCount)
 571         {
 572             ener_file_t    edrFile;
 573             gmx_enxnm_t   *energyTermsEdr = nullptr;
 574             int            numEnergyTermsEdr;
 575
 576             edrFile = open_enx(fileName, "r");
 577             do_enxnms(edrFile, &numEnergyTermsEdr, &energyTermsEdr);
 578             assert(energyTermsEdr);
 579
 580             // Check header
 581             TestReferenceChecker edrFileRef(checker_.checkCompound("File", "EnergyFile"));
 582             TestReferenceChecker energyTermsRef(edrFileRef.checkSequenceCompound("EnergyTerms", numEnergyTermsEdr));
 583             for (int i = 0; i < numEnergyTermsEdr; i++)
 584             {
 585                 TestReferenceChecker energyTermRef(energyTermsRef.checkCompound("EnergyTerm", nullptr));
 586                 energyTermRef.checkString(energyTermsEdr[i].name, "Name");
 587                 energyTermRef.checkString(energyTermsEdr[i].unit, "Units");
 588             }
 589
 590             // Check frames
 591             TestReferenceChecker framesRef(edrFileRef.checkSequenceCompound("Frames", frameCount));
 592             t_enxframe          *frameEdr;
 593             snew(frameEdr, 1);
 594             char                 buffer[22];
 595             for (int frameId = 0; frameId < frameCount; frameId++)
 596             {
 597                 bool bCont = do_enx(edrFile, frameEdr);
 598                 EXPECT_TRUE(bCont) << gmx::formatString("Cant read frame %d from .edr file.", frameId);
 599
 600                 TestReferenceChecker frameRef(framesRef.checkCompound("Frame", nullptr));
 601                 frameRef.checkReal(frameEdr->t, "Time");
 602                 frameRef.checkReal(frameEdr->dt, "Timestep");
 603                 frameRef.checkString(gmx_step_str(frameEdr->step, buffer), "Step");
 604                 frameRef.checkString(gmx_step_str(frameEdr->nsum, buffer), "NumSteps");
 605
 606                 EXPECT_EQ(frameEdr->nre, numEnergyTermsEdr) << gmx::formatString("Wrong number of energy terms in frame %d.", frameId);
 607                 TestReferenceChecker energyValuesRef(frameRef.checkSequenceCompound("EnergyTerms", numEnergyTermsEdr));
 608                 for (int i = 0; i < numEnergyTermsEdr; i++)
 609                 {
 610                     TestReferenceChecker energyValueRef(energyValuesRef.checkCompound("EnergyTerm", nullptr));
 611                     energyValueRef.checkString(energyTermsEdr[i].name, "Name");
 612                     energyValueRef.checkReal(frameEdr->ener[i].e, "Value");
 613                 }
 614             }
 615
 616             free_enxnms(numEnergyTermsEdr, energyTermsEdr);
 617             done_ener_file(edrFile);
 618
 619             free_enxframe(frameEdr);
 620             sfree(frameEdr);
 621         }
 622
 623 };
 624
 625 TEST_P(EnergyOutputTest, CheckOutput)
 626 {
 627     ASSERT_NE(log_, nullptr);
 628     // Binary output will be written to the temporary location
 629     energyFile_ = open_enx(edrFilename_.c_str(), "w");
 630     ASSERT_NE(energyFile_, nullptr);
 631
 632     EnergyOutputTestParameters parameters = GetParam();
 633     inputrec_.etc = parameters.temperatureCouplingScheme;
 634     inputrec_.epc = parameters.pressureCouplingScheme;
 635     inputrec_.eI  = parameters.integrator;
 636
 637     if (parameters.isBoxTriclinic)
 638     {
 639         inputrec_.ref_p[YY][XX]   = 1.0;
 640     }
 641
 642     energyOutput_.prepare(energyFile_, &mtop_, &inputrec_, nullptr, nullptr, parameters.isRerun);
 643
 644     // Add synthetic data for a single step
 645     double testValue = 10.0;
 646     for (int frame = 0; frame < parameters.numFrames; frame++)
 647     {
 648         setStepData(&testValue);
 649         energyOutput_.addDataAtEnergyStep(false, true, time_, tmass_, enerdata_.get(),
 650                                           &state_, nullptr, nullptr, box_,
 651                                           constraintsVirial_, forceVirial_, totalVirial_, pressure_,
 652                                           &ekindata_, muTotal_, constraints_.get());
 653
 654         energyOutput_.printStepToEnergyFile(energyFile_, true, false, false, log_,
 655                                             100*frame, time_, eprNORMAL,
 656                                             nullptr,  &mtop_.groups, &inputrec_.opts, nullptr);
 657         time_ += 1.0;
 658     }
 659
 660     energyOutput_.printStepToEnergyFile(energyFile_, true, false, false, log_,
 661                                         100, time_, eprAVER,
 662                                         nullptr,  &mtop_.groups, &inputrec_.opts, nullptr);
 663
 664     // We need to close the file before the contents are available.
 665     logFileGuard_.reset(nullptr);
 666
 667     done_ener_file(energyFile_);
 668
 669     // Set tolerance
 670     checker_.setDefaultTolerance(relativeToleranceAsFloatingPoint(testValue, 1.0e-5));
 671
 672     if (parameters.numFrames > 0)
 673     {
 674         // Test binary output
 675         checkEdrFile(edrFilename_.c_str(), parameters.numFrames);
 676     }
 677
 678     // Test printed values
 679     checker_.checkInteger(energyOutput_.numEnergyTerms(), "Number of Energy Terms");
 680     checker_.checkString(TextReader::readFileToString(logFilename_), "log");
 681 }
 682
 683 INSTANTIATE_TEST_CASE_P(WithParameters, EnergyOutputTest,
 684                             ::testing::ValuesIn(parametersSets));
 685
 686 }  // namespace
 687 }  // namespace test
 688 }  // namespace gmx