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