2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
37 * Tests for energy output to log and .edr files.
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.
47 * \author Mark Abraham <mark.j.abraham@gmail.com>
48 * \author Artem Zhmurov <zhmurov@gmail.com>
50 * \ingroup module_mdlib
54 #include "gromacs/mdlib/energyoutput.h"
58 #include <gtest/gtest.h>
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"
76 #include "testutils/refdata.h"
77 #include "testutils/setenv.h"
78 #include "testutils/testasserts.h"
79 #include "testutils/testfilemanager.h"
88 //! Wraps fclose to discard the return value to use it as a deleter with gmx::unique_cptr.
89 void fcloseWrapper(FILE* fp)
94 /*! \brief Test parameters space.
96 * The test will run on a set of combinations of this steucture parameters.
98 struct EnergyOutputTestParameters
100 //! Thermostat (enum)
101 TemperatureCoupling temperatureCouplingScheme;
103 PressureCoupling pressureCouplingScheme;
105 IntegrationAlgorithm integrator;
106 //! Number of saved energy frames (to test averages output).
108 //! If output should be initialized as a rerun.
110 //! Is box triclinic (off-diagonal elements will be printed).
114 /*! \brief Sets of parameters on which to run the tests.
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.
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 }
133 /*! \brief Test fixture to test energy output.
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.
139 class EnergyOutputTest : public ::testing::TestWithParam<EnergyOutputTestParameters>
141 int numTempCouplingGroups_ = 3;
142 real cosAccel_ = 1.0;
146 TestFileManager fileManager_;
147 //! Energy (.edr) file
148 ener_file_t energyFile_;
150 t_inputrec inputrec_;
157 //! Potential energy data
158 std::unique_ptr<gmx_enerdata_t> enerdata_;
159 //! Kinetic energy data (for temperatures output)
160 gmx_ekindata_t ekindata_;
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
177 //! Communication record
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
188 unique_cptr<FILE, fcloseWrapper> logFileGuard_;
190 TestReferenceData refData_;
191 //! Checker for reference data
192 TestReferenceChecker checker_;
195 ekindata_(numTempCouplingGroups_, cosAccel_, 1),
196 logFilename_(fileManager_.getTemporaryFilePath(".log")),
197 edrFilename_(fileManager_.getTemporaryFilePath(".edr")),
198 log_(std::fopen(logFilename_.c_str(), "w")),
200 checker_(refData_.rootChecker())
203 inputrec_.delta_t = 0.001;
205 // F_RF_EXCL will not be tested - group scheme is not supported any more
206 inputrec_.cutoff_scheme = CutoffScheme::Verlet;
208 inputrec_.coulombtype = CoulombInteractionType::Pme;
210 inputrec_.vdwtype = VanDerWaalsType::Pme;
212 // F_DVDL_COUL, F_DVDL_VDW, F_DVDL_BONDED, F_DVDL_RESTRAINT, F_DKDL and F_DVDL
213 inputrec_.efep = FreeEnergyPerturbationType::Yes;
214 inputrec_.fepvals->separate_dvdl[FreeEnergyPerturbationCouplingType::Coul] = true;
215 inputrec_.fepvals->separate_dvdl[FreeEnergyPerturbationCouplingType::Vdw] = true;
216 inputrec_.fepvals->separate_dvdl[FreeEnergyPerturbationCouplingType::Bonded] = true;
217 inputrec_.fepvals->separate_dvdl[FreeEnergyPerturbationCouplingType::Restraint] = true;
218 inputrec_.fepvals->separate_dvdl[FreeEnergyPerturbationCouplingType::Mass] = true;
219 inputrec_.fepvals->separate_dvdl[FreeEnergyPerturbationCouplingType::Coul] = true;
220 inputrec_.fepvals->separate_dvdl[FreeEnergyPerturbationCouplingType::Fep] = true;
222 // F_DISPCORR and F_PDISPCORR
223 inputrec_.eDispCorr = DispersionCorrectionType::Ener;
224 inputrec_.bRot = true;
227 inputrec_.ref_p[YY][XX] = 0.0;
228 inputrec_.ref_p[ZZ][XX] = 0.0;
229 inputrec_.ref_p[ZZ][YY] = 0.0;
232 inputrec_.ewald_geometry = EwaldGeometry::ThreeDC;
234 // To print constrain RMSD, constraints algorithm should be set to LINCS.
235 inputrec_.eConstrAlg = ConstraintAlgorithm::Lincs;
237 mtop_.bIntermolecularInteractions = false;
239 // Constructing molecular topology
240 gmx_moltype_t molType;
242 molType.atoms.nr = 2;
245 // This must be initialized so that Constraints object can be created below.
246 InteractionList interactionListConstr;
247 interactionListConstr.iatoms.resize(NRAL(F_CONSTR) + 1);
248 interactionListConstr.iatoms[0] = 0;
249 interactionListConstr.iatoms[1] = 0;
250 interactionListConstr.iatoms[2] = 1;
251 molType.ilist.at(F_CONSTR) = interactionListConstr;
253 InteractionList interactionListEmpty;
254 interactionListEmpty.iatoms.resize(0);
255 molType.ilist.at(F_CONSTRNC) = interactionListEmpty;
256 molType.ilist.at(F_SETTLE) = interactionListEmpty;
258 // F_LJ14 and F_COUL14
259 InteractionList interactionListLJ14;
260 interactionListLJ14.iatoms.resize(NRAL(F_LJ14) + 1);
261 molType.ilist.at(F_LJ14) = interactionListLJ14;
264 InteractionList interactionListLJC14Q;
265 interactionListLJC14Q.iatoms.resize(NRAL(F_LJC14_Q) + 1);
266 molType.ilist.at(F_LJC14_Q) = interactionListLJC14Q;
268 // TODO Do proper initialization for distance and orientation
269 // restraints and remove comments to enable their output
271 // InteractionList interactionListDISRES;
272 // interactionListDISRES.iatoms.resize(NRAL(F_DISRES) + 1);
273 // molType.ilist.at(F_DISRES) = interactionListDISRES;
276 // InteractionList interactionListORIRES;
277 // interactionListORIRES.iatoms.resize(NRAL(F_ORIRES) + 1);
278 // molType.ilist.at(F_ORIRES) = interactionListORIRES;
280 mtop_.moltype.push_back(molType);
282 gmx_molblock_t molBlock;
285 mtop_.molblock.push_back(molBlock);
287 // This is to keep constraints initialization happy
289 mtop_.ffparams.iparams.resize(F_NRE);
290 mtop_.ffparams.functype.resize(F_NRE);
291 mtop_.ffparams.iparams.at(F_CONSTR).constr.dA = 1.0;
292 mtop_.ffparams.iparams.at(F_CONSTR).constr.dB = 1.0;
293 mtop_.ffparams.iparams.at(F_CONSTRNC).constr.dA = 1.0;
294 mtop_.ffparams.iparams.at(F_CONSTRNC).constr.dB = 1.0;
296 // Groups for energy output, temperature coupling and acceleration
297 for (const auto& string : groupNameStrings_)
299 std::vector<char> cString(string.begin(), string.end());
300 // Need to add null termination
301 cString.push_back('\0');
302 groupNameCStrings_.emplace_back(cString);
303 groupNameHandles_.emplace_back(groupNameCStrings_.back().data());
305 for (auto& handle : groupNameHandles_)
307 mtop_.groups.groupNames.emplace_back(&handle);
310 mtop_.groups.groups[SimulationAtomGroupType::EnergyOutput].resize(numTempCouplingGroups_);
311 mtop_.groups.groups[SimulationAtomGroupType::EnergyOutput][0] = 0;
312 mtop_.groups.groups[SimulationAtomGroupType::EnergyOutput][1] = 1;
313 mtop_.groups.groups[SimulationAtomGroupType::EnergyOutput][2] = 2;
315 mtop_.groups.groups[SimulationAtomGroupType::TemperatureCoupling].resize(numTempCouplingGroups_);
316 mtop_.groups.groups[SimulationAtomGroupType::TemperatureCoupling][0] = 0;
317 mtop_.groups.groups[SimulationAtomGroupType::TemperatureCoupling][1] = 1;
318 mtop_.groups.groups[SimulationAtomGroupType::TemperatureCoupling][2] = 2;
320 // Nose-Hoover chains
321 inputrec_.bPrintNHChains = true;
322 inputrec_.opts.nhchainlength = 2;
323 state_.nosehoover_xi.resize(
324 mtop_.groups.groups[SimulationAtomGroupType::TemperatureCoupling].size()
325 * inputrec_.opts.nhchainlength);
326 state_.nosehoover_vxi.resize(
327 mtop_.groups.groups[SimulationAtomGroupType::TemperatureCoupling].size()
328 * inputrec_.opts.nhchainlength);
330 // This will be needed only with MTTK barostat
331 state_.nhpres_xi.resize(1 * inputrec_.opts.nhchainlength);
332 state_.nhpres_vxi.resize(1 * inputrec_.opts.nhchainlength);
335 enerdata_ = std::make_unique<gmx_enerdata_t>(
336 mtop_.groups.groups[SimulationAtomGroupType::EnergyOutput].size(), 0);
338 // Kinetic energy and related data
339 ekindata_.tcstat.resize(mtop_.groups.groups[SimulationAtomGroupType::TemperatureCoupling].size());
341 // This is needed so that the ebin space will be allocated
342 inputrec_.cos_accel = cosAccel_;
344 // Group options for annealing output
345 inputrec_.opts.ngtc = numTempCouplingGroups_;
346 snew(inputrec_.opts.ref_t, inputrec_.opts.ngtc);
347 snew(inputrec_.opts.annealing, inputrec_.opts.ngtc);
348 inputrec_.opts.annealing[0] = SimulatedAnnealing::No;
349 inputrec_.opts.annealing[1] = SimulatedAnnealing::Single;
350 inputrec_.opts.annealing[2] = SimulatedAnnealing::Periodic;
352 // This is to keep done_inputrec happy (otherwise sfree() segfaults)
353 snew(inputrec_.opts.anneal_time, inputrec_.opts.ngtc);
354 snew(inputrec_.opts.anneal_temp, inputrec_.opts.ngtc);
356 // Communication record (for Constraints constructor)
360 // Constraints object (to get constraints RMSD from)
361 // TODO EnergyOutput should not take Constraints object
362 // TODO This object will always return zero as RMSD value.
363 // It is more relevant to have non-zero value for testing.
364 constraints_ = makeConstraints(
365 mtop_, inputrec_, nullptr, false, nullptr, &cr_, false, nullptr, nullptr, nullptr, false, nullptr);
368 /*! \brief Helper function to generate synthetic data to output
370 * \param[in,out] testValue Base value fr energy data.
372 void setStepData(double* testValue)
375 time_ = (*testValue += 0.1);
376 tmass_ = (*testValue += 0.1);
378 enerdata_->term[F_LJ] = (*testValue += 0.1);
379 enerdata_->term[F_COUL_SR] = (*testValue += 0.1);
380 enerdata_->term[F_EPOT] = (*testValue += 0.1);
381 enerdata_->term[F_EKIN] = (*testValue += 0.1);
382 enerdata_->term[F_ETOT] = (*testValue += 0.1);
383 enerdata_->term[F_TEMP] = (*testValue += 0.1);
384 enerdata_->term[F_PRES] = (*testValue += 0.1);
386 enerdata_->term[F_BHAM] = (*testValue += 0.1);
387 enerdata_->term[F_EQM] = (*testValue += 0.1);
388 enerdata_->term[F_RF_EXCL] = (*testValue += 0.1);
389 enerdata_->term[F_COUL_RECIP] = (*testValue += 0.1);
390 enerdata_->term[F_LJ_RECIP] = (*testValue += 0.1);
391 enerdata_->term[F_LJ14] = (*testValue += 0.1);
392 enerdata_->term[F_COUL14] = (*testValue += 0.1);
393 enerdata_->term[F_LJC14_Q] = (*testValue += 0.1);
394 enerdata_->term[F_LJC_PAIRS_NB] = (*testValue += 0.1);
396 enerdata_->term[F_DVDL_COUL] = (*testValue += 0.1);
397 enerdata_->term[F_DVDL_VDW] = (*testValue += 0.1);
398 enerdata_->term[F_DVDL_BONDED] = (*testValue += 0.1);
399 enerdata_->term[F_DVDL_RESTRAINT] = (*testValue += 0.1);
400 enerdata_->term[F_DKDL] = (*testValue += 0.1);
401 enerdata_->term[F_DVDL] = (*testValue += 0.1);
403 enerdata_->term[F_DISPCORR] = (*testValue += 0.1);
404 enerdata_->term[F_PDISPCORR] = (*testValue += 0.1);
405 enerdata_->term[F_DISRESVIOL] = (*testValue += 0.1);
406 enerdata_->term[F_ORIRESDEV] = (*testValue += 0.1);
407 enerdata_->term[F_COM_PULL] = (*testValue += 0.1);
408 enerdata_->term[F_ECONSERVED] = (*testValue += 0.1);
411 for (int i = 0; i < enerdata_->grpp.nener; i++)
413 for (int k = 0; k < static_cast<int>(NonBondedEnergyTerms::Count); k++)
415 enerdata_->grpp.energyGroupPairTerms[k][i] = (*testValue += 0.1);
419 // Kinetic energy and related data
420 for (auto& tcstat : ekindata_.tcstat)
422 tcstat.T = (*testValue += 0.1);
423 tcstat.lambda = (*testValue += 0.1);
425 // Removing constant acceleration removed a total increment of 0.6
426 // To avoid unnecessary changes in reference data, we keep the increment
429 // This conditional is to check whether the ebin was allocated.
430 // Otherwise it will print cosacc data into the first bin.
431 if (inputrec_.cos_accel != 0)
433 ekindata_.cosacc.cos_accel = (*testValue += 0.1);
434 ekindata_.cosacc.vcos = (*testValue += 0.1);
437 state_.box[XX][XX] = (*testValue += 0.1);
438 state_.box[XX][YY] = (*testValue += 0.1);
439 state_.box[XX][ZZ] = (*testValue += 0.1);
440 state_.box[YY][XX] = (*testValue += 0.1);
441 state_.box[YY][YY] = (*testValue += 0.1);
442 state_.box[YY][ZZ] = (*testValue += 0.1);
443 state_.box[ZZ][XX] = (*testValue += 0.1);
444 state_.box[ZZ][YY] = (*testValue += 0.1);
445 state_.box[ZZ][ZZ] = (*testValue += 0.1);
447 box_[XX][XX] = (*testValue += 0.1);
448 box_[XX][YY] = (*testValue += 0.1);
449 box_[XX][ZZ] = (*testValue += 0.1);
450 box_[YY][XX] = (*testValue += 0.1);
451 box_[YY][YY] = (*testValue += 0.1);
452 box_[YY][ZZ] = (*testValue += 0.1);
453 box_[ZZ][XX] = (*testValue += 0.1);
454 box_[ZZ][YY] = (*testValue += 0.1);
455 box_[ZZ][ZZ] = (*testValue += 0.1);
457 // Removing GMX_CONSTRVIR removed a total increment of 1.8
458 // To avoid unnecessary changes in reference data, we keep the increment
461 totalVirial_[XX][XX] = (*testValue += 0.1);
462 totalVirial_[XX][YY] = (*testValue += 0.1);
463 totalVirial_[XX][ZZ] = (*testValue += 0.1);
464 totalVirial_[YY][XX] = (*testValue += 0.1);
465 totalVirial_[YY][YY] = (*testValue += 0.1);
466 totalVirial_[YY][ZZ] = (*testValue += 0.1);
467 totalVirial_[ZZ][XX] = (*testValue += 0.1);
468 totalVirial_[ZZ][YY] = (*testValue += 0.1);
469 totalVirial_[ZZ][ZZ] = (*testValue += 0.1);
471 pressure_[XX][XX] = (*testValue += 0.1);
472 pressure_[XX][YY] = (*testValue += 0.1);
473 pressure_[XX][ZZ] = (*testValue += 0.1);
474 pressure_[YY][XX] = (*testValue += 0.1);
475 pressure_[YY][YY] = (*testValue += 0.1);
476 pressure_[YY][ZZ] = (*testValue += 0.1);
477 pressure_[ZZ][XX] = (*testValue += 0.1);
478 pressure_[ZZ][YY] = (*testValue += 0.1);
479 pressure_[ZZ][ZZ] = (*testValue += 0.1);
481 muTotal_[XX] = (*testValue += 0.1);
482 muTotal_[YY] = (*testValue += 0.1);
483 muTotal_[ZZ] = (*testValue += 0.1);
485 state_.boxv[XX][XX] = (*testValue += 0.1);
486 state_.boxv[XX][YY] = (*testValue += 0.1);
487 state_.boxv[XX][ZZ] = (*testValue += 0.1);
488 state_.boxv[YY][XX] = (*testValue += 0.1);
489 state_.boxv[YY][YY] = (*testValue += 0.1);
490 state_.boxv[YY][ZZ] = (*testValue += 0.1);
491 state_.boxv[ZZ][XX] = (*testValue += 0.1);
492 state_.boxv[ZZ][YY] = (*testValue += 0.1);
493 state_.boxv[ZZ][ZZ] = (*testValue += 0.1);
495 for (int i = 0; i < inputrec_.opts.ngtc; i++)
497 inputrec_.opts.ref_t[i] = (*testValue += 0.1);
500 for (index k = 0; k < ssize(mtop_.groups.groups[SimulationAtomGroupType::TemperatureCoupling])
501 * inputrec_.opts.nhchainlength;
504 state_.nosehoover_xi[k] = (*testValue += 0.1);
505 state_.nosehoover_vxi[k] = (*testValue += 0.1);
507 for (int k = 0; k < inputrec_.opts.nhchainlength; k++)
509 state_.nhpres_xi[k] = (*testValue += 0.1);
510 state_.nhpres_vxi[k] = (*testValue += 0.1);
514 /*! \brief Check if the contents of the .edr file correspond to the reference data.
516 * The code below is based on the 'gmx dump' tool.
518 * \param[in] fileName Name of the file to check.
519 * \param[in] frameCount Number of frames to check.
521 void checkEdrFile(const char* fileName, int frameCount)
524 gmx_enxnm_t* energyTermsEdr = nullptr;
525 int numEnergyTermsEdr;
527 edrFile = open_enx(fileName, "r");
528 do_enxnms(edrFile, &numEnergyTermsEdr, &energyTermsEdr);
529 assert(energyTermsEdr);
532 TestReferenceChecker edrFileRef(checker_.checkCompound("File", "EnergyFile"));
533 TestReferenceChecker energyTermsRef(
534 edrFileRef.checkSequenceCompound("EnergyTerms", numEnergyTermsEdr));
535 for (int i = 0; i < numEnergyTermsEdr; i++)
537 TestReferenceChecker energyTermRef(energyTermsRef.checkCompound("EnergyTerm", nullptr));
538 energyTermRef.checkString(energyTermsEdr[i].name, "Name");
539 energyTermRef.checkString(energyTermsEdr[i].unit, "Units");
543 TestReferenceChecker framesRef(edrFileRef.checkSequenceCompound("Frames", frameCount));
544 t_enxframe* frameEdr;
547 for (int frameId = 0; frameId < frameCount; frameId++)
549 bool bCont = do_enx(edrFile, frameEdr);
550 EXPECT_TRUE(bCont) << gmx::formatString("Cant read frame %d from .edr file.", frameId);
552 TestReferenceChecker frameRef(framesRef.checkCompound("Frame", nullptr));
553 frameRef.checkReal(frameEdr->t, "Time");
554 frameRef.checkReal(frameEdr->dt, "Timestep");
555 frameRef.checkString(gmx_step_str(frameEdr->step, buffer), "Step");
556 frameRef.checkString(gmx_step_str(frameEdr->nsum, buffer), "NumSteps");
558 EXPECT_EQ(frameEdr->nre, numEnergyTermsEdr)
559 << gmx::formatString("Wrong number of energy terms in frame %d.", frameId);
560 TestReferenceChecker energyValuesRef(
561 frameRef.checkSequenceCompound("EnergyTerms", numEnergyTermsEdr));
562 for (int i = 0; i < numEnergyTermsEdr; i++)
564 TestReferenceChecker energyValueRef(energyValuesRef.checkCompound("EnergyTerm", nullptr));
565 energyValueRef.checkString(energyTermsEdr[i].name, "Name");
566 energyValueRef.checkReal(frameEdr->ener[i].e, "Value");
570 free_enxnms(numEnergyTermsEdr, energyTermsEdr);
571 done_ener_file(edrFile);
573 free_enxframe(frameEdr);
578 TEST_P(EnergyOutputTest, CheckOutput)
580 ASSERT_NE(log_, nullptr);
581 // Binary output will be written to the temporary location
582 energyFile_ = open_enx(edrFilename_.c_str(), "w");
583 ASSERT_NE(energyFile_, nullptr);
585 EnergyOutputTestParameters parameters = GetParam();
586 inputrec_.etc = parameters.temperatureCouplingScheme;
587 inputrec_.epc = parameters.pressureCouplingScheme;
588 inputrec_.eI = parameters.integrator;
590 if (parameters.isBoxTriclinic)
592 inputrec_.ref_p[YY][XX] = 1.0;
595 MDModulesNotifiers mdModulesNotifiers;
596 std::unique_ptr<EnergyOutput> energyOutput =
597 std::make_unique<EnergyOutput>(energyFile_,
603 StartingBehavior::NewSimulation,
607 // Add synthetic data for a single step
608 double testValue = 10.0;
609 for (int frame = 0; frame < parameters.numFrames; frame++)
611 setStepData(&testValue);
612 energyOutput->addDataAtEnergyStep(false,
620 PTCouplingArrays({ state_.boxv,
621 state_.nosehoover_xi,
622 state_.nosehoover_vxi,
624 state_.nhpres_vxi }),
632 energyOutput->printAnnealingTemperatures(log_, &mtop_.groups, &inputrec_.opts);
633 energyOutput->printStepToEnergyFile(
634 energyFile_, true, false, false, log_, 100 * frame, time_, nullptr, nullptr);
638 energyOutput->printAnnealingTemperatures(log_, &mtop_.groups, &inputrec_.opts);
639 energyOutput->printAverages(log_, &mtop_.groups);
641 // We need to close the file before the contents are available.
642 logFileGuard_.reset(nullptr);
644 done_ener_file(energyFile_);
647 checker_.setDefaultTolerance(relativeToleranceAsFloatingPoint(testValue, 1.0e-5));
649 if (parameters.numFrames > 0)
651 // Test binary output
652 checkEdrFile(edrFilename_.c_str(), parameters.numFrames);
655 // Test printed values
656 checker_.checkInteger(energyOutput->numEnergyTerms(), "Number of Energy Terms");
657 checker_.checkString(TextReader::readFileToString(logFilename_), "log");
660 INSTANTIATE_TEST_SUITE_P(WithParameters, EnergyOutputTest, ::testing::ValuesIn(parametersSets));