src/gromacs/mdlib/tests/energyoutput.cpp

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