1 /*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 2013,2014,2015,2016,2017 The GROMACS development team.
5 * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
6 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
7 * and including many others, as listed in the AUTHORS file in the
8 * top-level source directory and at http://www.gromacs.org.
9 *
10 * GROMACS is free software; you can redistribute it and/or
11 * modify it under the terms of the GNU Lesser General Public License
12 * as published by the Free Software Foundation; either version 2.1
13 * of the License, or (at your option) any later version.
14 *
15 * GROMACS is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 * Lesser General Public License for more details.
19 *
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with GROMACS; if not, see
22 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
23 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 *
25 * If you want to redistribute modifications to GROMACS, please
26 * consider that scientific software is very special. Version
27 * control is crucial - bugs must be traceable. We will be happy to
28 * consider code for inclusion in the official distribution, but
29 * derived work must not be called official GROMACS. Details are found
30 * in the README & COPYING files - if they are missing, get the
31 * official version at http://www.gromacs.org.
32 *
33 * To help us fund GROMACS development, we humbly ask that you cite
34 * the research papers on the package. Check out http://www.gromacs.org.
35 */
36 #include "gmxpre.h"
37
38 #include "mdoutf.h"
39
40 #include "config.h"
41
42 #include "gromacs/commandline/filenm.h"
43 #include "gromacs/domdec/collect.h"
44 #include "gromacs/domdec/domdec_struct.h"
45 #include "gromacs/fileio/checkpoint.h"
46 #include "gromacs/fileio/gmxfio.h"
47 #include "gromacs/fileio/tngio.h"
48 #include "gromacs/fileio/trrio.h"
49 #include "gromacs/fileio/xtcio.h"
50 #include "gromacs/fileio/xvgr.h"
51 #include "gromacs/math/vec.h"
52 #include "gromacs/mdlib/trajectory_writing.h"
53 #include "gromacs/mdrunutility/handlerestart.h"
54 #include "gromacs/mdrunutility/multisim.h"
55 #include "gromacs/mdtypes/awh_history.h"
56 #include "gromacs/mdtypes/commrec.h"
57 #include "gromacs/mdtypes/df_history.h"
58 #include "gromacs/mdtypes/edsamhistory.h"
59 #include "gromacs/mdtypes/energyhistory.h"
60 #include "gromacs/mdtypes/imdoutputprovider.h"
61 #include "gromacs/mdtypes/inputrec.h"
62 #include "gromacs/mdtypes/md_enums.h"
63 #include "gromacs/mdtypes/mdrunoptions.h"
64 #include "gromacs/mdtypes/observableshistory.h"
65 #include "gromacs/mdtypes/state.h"
66 #include "gromacs/mdtypes/swaphistory.h"
67 #include "gromacs/timing/wallcycle.h"
68 #include "gromacs/topology/topology.h"
69 #include "gromacs/utility/baseversion.h"
70 #include "gromacs/utility/fatalerror.h"
71 #include "gromacs/utility/pleasecite.h"
72 #include "gromacs/utility/programcontext.h"
73 #include "gromacs/utility/smalloc.h"
74 #include "gromacs/utility/sysinfo.h"
75
76 struct gmx_mdoutf
77 {
78 t_fileio* fp_trn;
79 t_fileio* fp_xtc;
80 gmx_tng_trajectory_t tng;
81 gmx_tng_trajectory_t tng_low_prec;
82 int x_compression_precision; /* only used by XTC output */
83 ener_file_t fp_ene;
84 const char* fn_cpt;
85 gmx_bool bKeepAndNumCPT;
86 int eIntegrator;
87 gmx_bool bExpanded;
88 int elamstats;
89 int simulation_part;
90 FILE* fp_dhdl;
91 int natoms_global;
92 int natoms_x_compressed;
93 const SimulationGroups* groups; /* for compressed position writing */
94 gmx_wallcycle_t wcycle;
95 rvec* f_global;
96 gmx::IMDOutputProvider* outputProvider;
97 const gmx::MdModulesNotifier* mdModulesNotifier;
98 bool simulationsShareState;
99 MPI_Comm mastersComm;
100 };
101
102
init_mdoutf(FILE * fplog,int nfile,const t_filenm fnm[],const gmx::MdrunOptions & mdrunOptions,const t_commrec * cr,gmx::IMDOutputProvider * outputProvider,const gmx::MdModulesNotifier & mdModulesNotifier,const t_inputrec * ir,const gmx_mtop_t * top_global,const gmx_output_env_t * oenv,gmx_wallcycle_t wcycle,const gmx::StartingBehavior startingBehavior,bool simulationsShareState,const gmx_multisim_t * ms)103 gmx_mdoutf_t init_mdoutf(FILE* fplog,
104 int nfile,
105 const t_filenm fnm[],
106 const gmx::MdrunOptions& mdrunOptions,
107 const t_commrec* cr,
108 gmx::IMDOutputProvider* outputProvider,
109 const gmx::MdModulesNotifier& mdModulesNotifier,
110 const t_inputrec* ir,
111 const gmx_mtop_t* top_global,
112 const gmx_output_env_t* oenv,
113 gmx_wallcycle_t wcycle,
114 const gmx::StartingBehavior startingBehavior,
115 bool simulationsShareState,
116 const gmx_multisim_t* ms)
117 {
118 gmx_mdoutf_t of;
119 const char * appendMode = "a+", *writeMode = "w+", *filemode;
120 gmx_bool bCiteTng = FALSE;
121 int i;
122 bool restartWithAppending = (startingBehavior == gmx::StartingBehavior::RestartWithAppending);
123
124 snew(of, 1);
125
126 of->fp_trn = nullptr;
127 of->fp_ene = nullptr;
128 of->fp_xtc = nullptr;
129 of->tng = nullptr;
130 of->tng_low_prec = nullptr;
131 of->fp_dhdl = nullptr;
132
133 of->eIntegrator = ir->eI;
134 of->bExpanded = ir->bExpanded;
135 of->elamstats = ir->expandedvals->elamstats;
136 of->simulation_part = ir->simulation_part;
137 of->x_compression_precision = static_cast<int>(ir->x_compression_precision);
138 of->wcycle = wcycle;
139 of->f_global = nullptr;
140 of->outputProvider = outputProvider;
141
142 GMX_RELEASE_ASSERT(!simulationsShareState || ms != nullptr,
143 "Need valid multisim object when simulations share state");
144 of->simulationsShareState = simulationsShareState;
145 if (of->simulationsShareState)
146 {
147 of->mastersComm = ms->mastersComm_;
148 }
149
150 if (MASTER(cr))
151 {
152 of->bKeepAndNumCPT = mdrunOptions.checkpointOptions.keepAndNumberCheckpointFiles;
153
154 filemode = restartWithAppending ? appendMode : writeMode;
155
156 if (EI_DYNAMICS(ir->eI) && ir->nstxout_compressed > 0)
157 {
158 const char* filename;
159 filename = ftp2fn(efCOMPRESSED, nfile, fnm);
160 switch (fn2ftp(filename))
161 {
162 case efXTC: of->fp_xtc = open_xtc(filename, filemode); break;
163 case efTNG:
164 gmx_tng_open(filename, filemode[0], &of->tng_low_prec);
165 if (filemode[0] == 'w')
166 {
167 gmx_tng_prepare_low_prec_writing(of->tng_low_prec, top_global, ir);
168 }
169 bCiteTng = TRUE;
170 break;
171 default: gmx_incons("Invalid reduced precision file format");
172 }
173 }
174 if ((EI_DYNAMICS(ir->eI) || EI_ENERGY_MINIMIZATION(ir->eI))
175 && (!GMX_FAHCORE
176 && !(EI_DYNAMICS(ir->eI) && ir->nstxout == 0 && ir->nstvout == 0 && ir->nstfout == 0)))
177 {
178 const char* filename;
179 filename = ftp2fn(efTRN, nfile, fnm);
180 switch (fn2ftp(filename))
181 {
182 case efTRR:
183 case efTRN:
184 /* If there is no uncompressed coordinate output and
185 there is compressed TNG output write forces
186 and/or velocities to the TNG file instead. */
187 if (ir->nstxout != 0 || ir->nstxout_compressed == 0 || !of->tng_low_prec)
188 {
189 of->fp_trn = gmx_trr_open(filename, filemode);
190 }
191 break;
192 case efTNG:
193 gmx_tng_open(filename, filemode[0], &of->tng);
194 if (filemode[0] == 'w')
195 {
196 gmx_tng_prepare_md_writing(of->tng, top_global, ir);
197 }
198 bCiteTng = TRUE;
199 break;
200 default: gmx_incons("Invalid full precision file format");
201 }
202 }
203 if (EI_DYNAMICS(ir->eI) || EI_ENERGY_MINIMIZATION(ir->eI))
204 {
205 of->fp_ene = open_enx(ftp2fn(efEDR, nfile, fnm), filemode);
206 }
207 of->fn_cpt = opt2fn("-cpo", nfile, fnm);
208
209 if ((ir->efep != efepNO || ir->bSimTemp) && ir->fepvals->nstdhdl > 0
210 && (ir->fepvals->separate_dhdl_file == esepdhdlfileYES) && EI_DYNAMICS(ir->eI))
211 {
212 if (restartWithAppending)
213 {
214 of->fp_dhdl = gmx_fio_fopen(opt2fn("-dhdl", nfile, fnm), filemode);
215 }
216 else
217 {
218 of->fp_dhdl = open_dhdl(opt2fn("-dhdl", nfile, fnm), ir, oenv);
219 }
220 }
221
222 outputProvider->initOutput(fplog, nfile, fnm, restartWithAppending, oenv);
223 of->mdModulesNotifier = &mdModulesNotifier;
224
225 /* Set up atom counts so they can be passed to actual
226 trajectory-writing routines later. Also, XTC writing needs
227 to know what (and how many) atoms might be in the XTC
228 groups, and how to look up later which ones they are. */
229 of->natoms_global = top_global->natoms;
230 of->groups = &top_global->groups;
231 of->natoms_x_compressed = 0;
232 for (i = 0; (i < top_global->natoms); i++)
233 {
234 if (getGroupType(*of->groups, SimulationAtomGroupType::CompressedPositionOutput, i) == 0)
235 {
236 of->natoms_x_compressed++;
237 }
238 }
239
240 if (ir->nstfout && DOMAINDECOMP(cr))
241 {
242 snew(of->f_global, top_global->natoms);
243 }
244 }
245
246 if (bCiteTng)
247 {
248 please_cite(fplog, "Lundborg2014");
249 }
250
251 return of;
252 }
253
mdoutf_get_fp_ene(gmx_mdoutf_t of)254 ener_file_t mdoutf_get_fp_ene(gmx_mdoutf_t of)
255 {
256 return of->fp_ene;
257 }
258
mdoutf_get_fp_dhdl(gmx_mdoutf_t of)259 FILE* mdoutf_get_fp_dhdl(gmx_mdoutf_t of)
260 {
261 return of->fp_dhdl;
262 }
263
mdoutf_get_wcycle(gmx_mdoutf_t of)264 gmx_wallcycle_t mdoutf_get_wcycle(gmx_mdoutf_t of)
265 {
266 return of->wcycle;
267 }
268
mpiBarrierBeforeRename(const bool applyMpiBarrierBeforeRename,MPI_Comm mpiBarrierCommunicator)269 static void mpiBarrierBeforeRename(const bool applyMpiBarrierBeforeRename, MPI_Comm mpiBarrierCommunicator)
270 {
271 if (applyMpiBarrierBeforeRename)
272 {
273 #if GMX_MPI
274 MPI_Barrier(mpiBarrierCommunicator);
275 #else
276 GMX_RELEASE_ASSERT(false, "Should not request a barrier without MPI");
277 GMX_UNUSED_VALUE(mpiBarrierCommunicator);
278 #endif
279 }
280 }
281 /*! \brief Write a checkpoint to the filename
282 *
283 * Appends the _step<step>.cpt with bNumberAndKeep, otherwise moves
284 * the previous checkpoint filename with suffix _prev.cpt.
285 */
write_checkpoint(const char * fn,gmx_bool bNumberAndKeep,FILE * fplog,const t_commrec * cr,ivec domdecCells,int nppnodes,int eIntegrator,int simulation_part,gmx_bool bExpanded,int elamstats,int64_t step,double t,t_state * state,ObservablesHistory * observablesHistory,const gmx::MdModulesNotifier & mdModulesNotifier,gmx::WriteCheckpointDataHolder * modularSimulatorCheckpointData,bool applyMpiBarrierBeforeRename,MPI_Comm mpiBarrierCommunicator)286 static void write_checkpoint(const char* fn,
287 gmx_bool bNumberAndKeep,
288 FILE* fplog,
289 const t_commrec* cr,
290 ivec domdecCells,
291 int nppnodes,
292 int eIntegrator,
293 int simulation_part,
294 gmx_bool bExpanded,
295 int elamstats,
296 int64_t step,
297 double t,
298 t_state* state,
299 ObservablesHistory* observablesHistory,
300 const gmx::MdModulesNotifier& mdModulesNotifier,
301 gmx::WriteCheckpointDataHolder* modularSimulatorCheckpointData,
302 bool applyMpiBarrierBeforeRename,
303 MPI_Comm mpiBarrierCommunicator)
304 {
305 t_fileio* fp;
306 char* fntemp; /* the temporary checkpoint file name */
307 int npmenodes;
308 char buf[1024], suffix[5 + STEPSTRSIZE], sbuf[STEPSTRSIZE];
309 t_fileio* ret;
310
311 if (DOMAINDECOMP(cr))
312 {
313 npmenodes = cr->npmenodes;
314 }
315 else
316 {
317 npmenodes = 0;
318 }
319
320 #if !GMX_NO_RENAME
321 /* make the new temporary filename */
322 snew(fntemp, std::strlen(fn) + 5 + STEPSTRSIZE);
323 std::strcpy(fntemp, fn);
324 fntemp[std::strlen(fn) - std::strlen(ftp2ext(fn2ftp(fn))) - 1] = '\0';
325 sprintf(suffix, "_%s%s", "step", gmx_step_str(step, sbuf));
326 std::strcat(fntemp, suffix);
327 std::strcat(fntemp, fn + std::strlen(fn) - std::strlen(ftp2ext(fn2ftp(fn))) - 1);
328 #else
329 /* if we can't rename, we just overwrite the cpt file.
330 * dangerous if interrupted.
331 */
332 snew(fntemp, std::strlen(fn));
333 std::strcpy(fntemp, fn);
334 #endif
335 std::string timebuf = gmx_format_current_time();
336
337 if (fplog)
338 {
339 fprintf(fplog, "Writing checkpoint, step %s at %s\n\n", gmx_step_str(step, buf), timebuf.c_str());
340 }
341
342 /* Get offsets for open files */
343 auto outputfiles = gmx_fio_get_output_file_positions();
344
345 fp = gmx_fio_open(fntemp, "w");
346
347 /* We can check many more things now (CPU, acceleration, etc), but
348 * it is highly unlikely to have two separate builds with exactly
349 * the same version, user, time, and build host!
350 */
351
352 int nlambda = (state->dfhist ? state->dfhist->nlambda : 0);
353
354 edsamhistory_t* edsamhist = observablesHistory->edsamHistory.get();
355 int nED = (edsamhist ? edsamhist->nED : 0);
356
357 swaphistory_t* swaphist = observablesHistory->swapHistory.get();
358 int eSwapCoords = (swaphist ? swaphist->eSwapCoords : eswapNO);
359
360 CheckpointHeaderContents headerContents = { 0,
361 { 0 },
362 { 0 },
363 { 0 },
364 { 0 },
365 GMX_DOUBLE,
366 { 0 },
367 { 0 },
368 eIntegrator,
369 simulation_part,
370 step,
371 t,
372 nppnodes,
373 { 0 },
374 npmenodes,
375 state->natoms,
376 state->ngtc,
377 state->nnhpres,
378 state->nhchainlength,
379 nlambda,
380 state->flags,
381 0,
382 0,
383 0,
384 0,
385 0,
386 nED,
387 eSwapCoords,
388 false };
389 std::strcpy(headerContents.version, gmx_version());
390 std::strcpy(headerContents.fprog, gmx::getProgramContext().fullBinaryPath());
391 std::strcpy(headerContents.ftime, timebuf.c_str());
392 if (DOMAINDECOMP(cr))
393 {
394 copy_ivec(domdecCells, headerContents.dd_nc);
395 }
396
397 write_checkpoint_data(fp, headerContents, bExpanded, elamstats, state, observablesHistory,
398 mdModulesNotifier, &outputfiles, modularSimulatorCheckpointData);
399
400 /* we really, REALLY, want to make sure to physically write the checkpoint,
401 and all the files it depends on, out to disk. Because we've
402 opened the checkpoint with gmx_fio_open(), it's in our list
403 of open files. */
404 ret = gmx_fio_all_output_fsync();
405
406 if (ret)
407 {
408 char buf[STRLEN];
409 sprintf(buf, "Cannot fsync '%s'; maybe you are out of disk space?", gmx_fio_getname(ret));
410
411 if (getenv(GMX_IGNORE_FSYNC_FAILURE_ENV) == nullptr)
412 {
413 gmx_file(buf);
414 }
415 else
416 {
417 gmx_warning("%s", buf);
418 }
419 }
420
421 if (gmx_fio_close(fp) != 0)
422 {
423 gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");
424 }
425
426 /* we don't move the checkpoint if the user specified they didn't want it,
427 or if the fsyncs failed */
428 #if !GMX_NO_RENAME
429 if (!bNumberAndKeep && !ret)
430 {
431 if (gmx_fexist(fn))
432 {
433 /* Rename the previous checkpoint file */
434 mpiBarrierBeforeRename(applyMpiBarrierBeforeRename, mpiBarrierCommunicator);
435
436 std::strcpy(buf, fn);
437 buf[std::strlen(fn) - std::strlen(ftp2ext(fn2ftp(fn))) - 1] = '\0';
438 std::strcat(buf, "_prev");
439 std::strcat(buf, fn + std::strlen(fn) - std::strlen(ftp2ext(fn2ftp(fn))) - 1);
440 if (!GMX_FAHCORE)
441 {
442 /* we copy here so that if something goes wrong between now and
443 * the rename below, there's always a state.cpt.
444 * If renames are atomic (such as in POSIX systems),
445 * this copying should be unneccesary.
446 */
447 gmx_file_copy(fn, buf, FALSE);
448 /* We don't really care if this fails:
449 * there's already a new checkpoint.
450 */
451 }
452 else
453 {
454 gmx_file_rename(fn, buf);
455 }
456 }
457
458 /* Rename the checkpoint file from the temporary to the final name */
459 mpiBarrierBeforeRename(applyMpiBarrierBeforeRename, mpiBarrierCommunicator);
460
461 if (gmx_file_rename(fntemp, fn) != 0)
462 {
463 gmx_file("Cannot rename checkpoint file; maybe you are out of disk space?");
464 }
465 }
466 #endif /* GMX_NO_RENAME */
467
468 sfree(fntemp);
469
470 #if GMX_FAHCORE
471 /* Always FAH checkpoint immediately after a GROMACS checkpoint.
472 *
473 * Note that it is critical that we save a FAH checkpoint directly
474 * after writing a GROMACS checkpoint. If the program dies, either
475 * by the machine powering off suddenly or the process being,
476 * killed, FAH can recover files that have only appended data by
477 * truncating them to the last recorded length. The GROMACS
478 * checkpoint does not just append data, it is fully rewritten each
479 * time so a crash between moving the new Gromacs checkpoint file in
480 * to place and writing a FAH checkpoint is not recoverable. Thus
481 * the time between these operations must be kept as short as
482 * possible.
483 */
484 fcCheckpoint();
485 #endif /* end GMX_FAHCORE block */
486 }
487
mdoutf_write_checkpoint(gmx_mdoutf_t of,FILE * fplog,const t_commrec * cr,int64_t step,double t,t_state * state_global,ObservablesHistory * observablesHistory,gmx::WriteCheckpointDataHolder * modularSimulatorCheckpointData)488 void mdoutf_write_checkpoint(gmx_mdoutf_t of,
489 FILE* fplog,
490 const t_commrec* cr,
491 int64_t step,
492 double t,
493 t_state* state_global,
494 ObservablesHistory* observablesHistory,
495 gmx::WriteCheckpointDataHolder* modularSimulatorCheckpointData)
496 {
497 fflush_tng(of->tng);
498 fflush_tng(of->tng_low_prec);
499 /* Write the checkpoint file.
500 * When simulations share the state, an MPI barrier is applied before
501 * renaming old and new checkpoint files to minimize the risk of
502 * checkpoint files getting out of sync.
503 */
504 ivec one_ivec = { 1, 1, 1 };
505 write_checkpoint(of->fn_cpt, of->bKeepAndNumCPT, fplog, cr,
506 DOMAINDECOMP(cr) ? cr->dd->numCells : one_ivec,
507 DOMAINDECOMP(cr) ? cr->dd->nnodes : cr->nnodes, of->eIntegrator,
508 of->simulation_part, of->bExpanded, of->elamstats, step, t, state_global,
509 observablesHistory, *(of->mdModulesNotifier), modularSimulatorCheckpointData,
510 of->simulationsShareState, of->mastersComm);
511 }
512
mdoutf_write_to_trajectory_files(FILE * fplog,const t_commrec * cr,gmx_mdoutf_t of,int mdof_flags,int natoms,int64_t step,double t,t_state * state_local,t_state * state_global,ObservablesHistory * observablesHistory,gmx::ArrayRef<const gmx::RVec> f_local,gmx::WriteCheckpointDataHolder * modularSimulatorCheckpointData)513 void mdoutf_write_to_trajectory_files(FILE* fplog,
514 const t_commrec* cr,
515 gmx_mdoutf_t of,
516 int mdof_flags,
517 int natoms,
518 int64_t step,
519 double t,
520 t_state* state_local,
521 t_state* state_global,
522 ObservablesHistory* observablesHistory,
523 gmx::ArrayRef<const gmx::RVec> f_local,
524 gmx::WriteCheckpointDataHolder* modularSimulatorCheckpointData)
525 {
526 const rvec* f_global;
527
528 if (DOMAINDECOMP(cr))
529 {
530 if (mdof_flags & MDOF_CPT)
531 {
532 dd_collect_state(cr->dd, state_local, state_global);
533 }
534 else
535 {
536 if (mdof_flags & (MDOF_X | MDOF_X_COMPRESSED))
537 {
538 auto globalXRef = MASTER(cr) ? state_global->x : gmx::ArrayRef<gmx::RVec>();
539 dd_collect_vec(cr->dd, state_local->ddp_count, state_local->ddp_count_cg_gl,
540 state_local->cg_gl, state_local->x, globalXRef);
541 }
542 if (mdof_flags & MDOF_V)
543 {
544 auto globalVRef = MASTER(cr) ? state_global->v : gmx::ArrayRef<gmx::RVec>();
545 dd_collect_vec(cr->dd, state_local->ddp_count, state_local->ddp_count_cg_gl,
546 state_local->cg_gl, state_local->v, globalVRef);
547 }
548 }
549 f_global = of->f_global;
550 if (mdof_flags & MDOF_F)
551 {
552 dd_collect_vec(
553 cr->dd, state_local->ddp_count, state_local->ddp_count_cg_gl, state_local->cg_gl, f_local,
554 gmx::arrayRefFromArray(reinterpret_cast<gmx::RVec*>(of->f_global), f_local.size()));
555 }
556 }
557 else
558 {
559 /* We have the whole state locally: copy the local state pointer */
560 state_global = state_local;
561
562 f_global = as_rvec_array(f_local.data());
563 }
564
565 if (MASTER(cr))
566 {
567 if (mdof_flags & MDOF_CPT)
568 {
569 mdoutf_write_checkpoint(of, fplog, cr, step, t, state_global, observablesHistory,
570 modularSimulatorCheckpointData);
571 }
572
573 if (mdof_flags & (MDOF_X | MDOF_V | MDOF_F))
574 {
575 const rvec* x = (mdof_flags & MDOF_X) ? state_global->x.rvec_array() : nullptr;
576 const rvec* v = (mdof_flags & MDOF_V) ? state_global->v.rvec_array() : nullptr;
577 const rvec* f = (mdof_flags & MDOF_F) ? f_global : nullptr;
578
579 if (of->fp_trn)
580 {
581 gmx_trr_write_frame(of->fp_trn, step, t, state_local->lambda[efptFEP],
582 state_local->box, natoms, x, v, f);
583 if (gmx_fio_flush(of->fp_trn) != 0)
584 {
585 gmx_file("Cannot write trajectory; maybe you are out of disk space?");
586 }
587 }
588
589 /* If a TNG file is open for uncompressed coordinate output also write
590 velocities and forces to it. */
591 else if (of->tng)
592 {
593 gmx_fwrite_tng(of->tng, FALSE, step, t, state_local->lambda[efptFEP],
594 state_local->box, natoms, x, v, f);
595 }
596 /* If only a TNG file is open for compressed coordinate output (no uncompressed
597 coordinate output) also write forces and velocities to it. */
598 else if (of->tng_low_prec)
599 {
600 gmx_fwrite_tng(of->tng_low_prec, FALSE, step, t, state_local->lambda[efptFEP],
601 state_local->box, natoms, x, v, f);
602 }
603 }
604 if (mdof_flags & MDOF_X_COMPRESSED)
605 {
606 rvec* xxtc = nullptr;
607
608 if (of->natoms_x_compressed == of->natoms_global)
609 {
610 /* We are writing the positions of all of the atoms to
611 the compressed output */
612 xxtc = state_global->x.rvec_array();
613 }
614 else
615 {
616 /* We are writing the positions of only a subset of
617 the atoms to the compressed output, so we have to
618 make a copy of the subset of coordinates. */
619 int i, j;
620
621 snew(xxtc, of->natoms_x_compressed);
622 auto x = makeArrayRef(state_global->x);
623 for (i = 0, j = 0; (i < of->natoms_global); i++)
624 {
625 if (getGroupType(*of->groups, SimulationAtomGroupType::CompressedPositionOutput, i) == 0)
626 {
627 copy_rvec(x[i], xxtc[j++]);
628 }
629 }
630 }
631 if (write_xtc(of->fp_xtc, of->natoms_x_compressed, step, t, state_local->box, xxtc,
632 of->x_compression_precision)
633 == 0)
634 {
635 gmx_fatal(FARGS,
636 "XTC error. This indicates you are out of disk space, or a "
637 "simulation with major instabilities resulting in coordinates "
638 "that are NaN or too large to be represented in the XTC format.\n");
639 }
640 gmx_fwrite_tng(of->tng_low_prec, TRUE, step, t, state_local->lambda[efptFEP],
641 state_local->box, of->natoms_x_compressed, xxtc, nullptr, nullptr);
642 if (of->natoms_x_compressed != of->natoms_global)
643 {
644 sfree(xxtc);
645 }
646 }
647 if (mdof_flags & (MDOF_BOX | MDOF_LAMBDA) && !(mdof_flags & (MDOF_X | MDOF_V | MDOF_F)))
648 {
649 if (of->tng)
650 {
651 real lambda = -1;
652 rvec* box = nullptr;
653 if (mdof_flags & MDOF_BOX)
654 {
655 box = state_local->box;
656 }
657 if (mdof_flags & MDOF_LAMBDA)
658 {
659 lambda = state_local->lambda[efptFEP];
660 }
661 gmx_fwrite_tng(of->tng, FALSE, step, t, lambda, box, natoms, nullptr, nullptr, nullptr);
662 }
663 }
664 if (mdof_flags & (MDOF_BOX_COMPRESSED | MDOF_LAMBDA_COMPRESSED)
665 && !(mdof_flags & (MDOF_X_COMPRESSED)))
666 {
667 if (of->tng_low_prec)
668 {
669 real lambda = -1;
670 rvec* box = nullptr;
671 if (mdof_flags & MDOF_BOX_COMPRESSED)
672 {
673 box = state_local->box;
674 }
675 if (mdof_flags & MDOF_LAMBDA_COMPRESSED)
676 {
677 lambda = state_local->lambda[efptFEP];
678 }
679 gmx_fwrite_tng(of->tng_low_prec, FALSE, step, t, lambda, box, natoms, nullptr,
680 nullptr, nullptr);
681 }
682 }
683
684 #if GMX_FAHCORE
685 /* Write a FAH checkpoint after writing any other data. We may end up
686 checkpointing twice but it's fast so it's ok. */
687 if ((mdof_flags & ~MDOF_CPT))
688 {
689 fcCheckpoint();
690 }
691 #endif
692 }
693 }
694
mdoutf_tng_close(gmx_mdoutf_t of)695 void mdoutf_tng_close(gmx_mdoutf_t of)
696 {
697 if (of->tng || of->tng_low_prec)
698 {
699 wallcycle_start(of->wcycle, ewcTRAJ);
700 gmx_tng_close(&of->tng);
701 gmx_tng_close(&of->tng_low_prec);
702 wallcycle_stop(of->wcycle, ewcTRAJ);
703 }
704 }
705
done_mdoutf(gmx_mdoutf_t of)706 void done_mdoutf(gmx_mdoutf_t of)
707 {
708 if (of->fp_ene != nullptr)
709 {
710 done_ener_file(of->fp_ene);
711 }
712 if (of->fp_xtc)
713 {
714 close_xtc(of->fp_xtc);
715 }
716 if (of->fp_trn)
717 {
718 gmx_trr_close(of->fp_trn);
719 }
720 if (of->fp_dhdl != nullptr)
721 {
722 gmx_fio_fclose(of->fp_dhdl);
723 }
724 of->outputProvider->finishOutput();
725 if (of->f_global != nullptr)
726 {
727 sfree(of->f_global);
728 }
729
730 gmx_tng_close(&of->tng);
731 gmx_tng_close(&of->tng_low_prec);
732
733 sfree(of);
734 }
735
mdoutf_get_tng_box_output_interval(gmx_mdoutf_t of)736 int mdoutf_get_tng_box_output_interval(gmx_mdoutf_t of)
737 {
738 if (of->tng)
739 {
740 return gmx_tng_get_box_output_interval(of->tng);
741 }
742 return 0;
743 }
744
mdoutf_get_tng_lambda_output_interval(gmx_mdoutf_t of)745 int mdoutf_get_tng_lambda_output_interval(gmx_mdoutf_t of)
746 {
747 if (of->tng)
748 {
749 return gmx_tng_get_lambda_output_interval(of->tng);
750 }
751 return 0;
752 }
753
mdoutf_get_tng_compressed_box_output_interval(gmx_mdoutf_t of)754 int mdoutf_get_tng_compressed_box_output_interval(gmx_mdoutf_t of)
755 {
756 if (of->tng_low_prec)
757 {
758 return gmx_tng_get_box_output_interval(of->tng_low_prec);
759 }
760 return 0;
761 }
762
mdoutf_get_tng_compressed_lambda_output_interval(gmx_mdoutf_t of)763 int mdoutf_get_tng_compressed_lambda_output_interval(gmx_mdoutf_t of)
764 {
765 if (of->tng_low_prec)
766 {
767 return gmx_tng_get_lambda_output_interval(of->tng_low_prec);
768 }
769 return 0;
770 }
771