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