2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2013,2014,2015,2016,2017 The GROMACS development team.
5 * Copyright (c) 2018,2019,2020, 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.
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.
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.
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.
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.
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.
40 #include "gromacs/commandline/filenm.h"
41 #include "gromacs/domdec/collect.h"
42 #include "gromacs/domdec/domdec_struct.h"
43 #include "gromacs/fileio/checkpoint.h"
44 #include "gromacs/fileio/gmxfio.h"
45 #include "gromacs/fileio/tngio.h"
46 #include "gromacs/fileio/trrio.h"
47 #include "gromacs/fileio/xtcio.h"
48 #include "gromacs/fileio/xvgr.h"
49 #include "gromacs/math/vec.h"
50 #include "gromacs/mdlib/trajectory_writing.h"
51 #include "gromacs/mdrunutility/handlerestart.h"
52 #include "gromacs/mdrunutility/multisim.h"
53 #include "gromacs/mdtypes/commrec.h"
54 #include "gromacs/mdtypes/imdoutputprovider.h"
55 #include "gromacs/mdtypes/inputrec.h"
56 #include "gromacs/mdtypes/md_enums.h"
57 #include "gromacs/mdtypes/mdrunoptions.h"
58 #include "gromacs/mdtypes/state.h"
59 #include "gromacs/timing/wallcycle.h"
60 #include "gromacs/topology/topology.h"
61 #include "gromacs/utility/fatalerror.h"
62 #include "gromacs/utility/pleasecite.h"
63 #include "gromacs/utility/smalloc.h"
69 gmx_tng_trajectory_t tng;
70 gmx_tng_trajectory_t tng_low_prec;
71 int x_compression_precision; /* only used by XTC output */
74 gmx_bool bKeepAndNumCPT;
81 int natoms_x_compressed;
82 const SimulationGroups* groups; /* for compressed position writing */
83 gmx_wallcycle_t wcycle;
85 gmx::IMDOutputProvider* outputProvider;
86 const gmx::MdModulesNotifier* mdModulesNotifier;
87 bool simulationsShareState;
88 MPI_Comm mpiCommMasters;
92 gmx_mdoutf_t init_mdoutf(FILE* fplog,
95 const gmx::MdrunOptions& mdrunOptions,
97 gmx::IMDOutputProvider* outputProvider,
98 const gmx::MdModulesNotifier& mdModulesNotifier,
100 const gmx_mtop_t* top_global,
101 const gmx_output_env_t* oenv,
102 gmx_wallcycle_t wcycle,
103 const gmx::StartingBehavior startingBehavior,
104 bool simulationsShareState,
105 const gmx_multisim_t* ms)
108 const char * appendMode = "a+", *writeMode = "w+", *filemode;
109 gmx_bool bCiteTng = FALSE;
111 bool restartWithAppending = (startingBehavior == gmx::StartingBehavior::RestartWithAppending);
115 of->fp_trn = nullptr;
116 of->fp_ene = nullptr;
117 of->fp_xtc = nullptr;
119 of->tng_low_prec = nullptr;
120 of->fp_dhdl = nullptr;
122 of->eIntegrator = ir->eI;
123 of->bExpanded = ir->bExpanded;
124 of->elamstats = ir->expandedvals->elamstats;
125 of->simulation_part = ir->simulation_part;
126 of->x_compression_precision = static_cast<int>(ir->x_compression_precision);
128 of->f_global = nullptr;
129 of->outputProvider = outputProvider;
131 GMX_RELEASE_ASSERT(!simulationsShareState || ms != nullptr,
132 "Need valid multisim object when simulations share state");
133 of->simulationsShareState = simulationsShareState;
134 if (of->simulationsShareState)
136 of->mpiCommMasters = ms->mpi_comm_masters;
141 of->bKeepAndNumCPT = mdrunOptions.checkpointOptions.keepAndNumberCheckpointFiles;
143 filemode = restartWithAppending ? appendMode : writeMode;
145 if (EI_DYNAMICS(ir->eI) && ir->nstxout_compressed > 0)
147 const char* filename;
148 filename = ftp2fn(efCOMPRESSED, nfile, fnm);
149 switch (fn2ftp(filename))
151 case efXTC: of->fp_xtc = open_xtc(filename, filemode); break;
153 gmx_tng_open(filename, filemode[0], &of->tng_low_prec);
154 if (filemode[0] == 'w')
156 gmx_tng_prepare_low_prec_writing(of->tng_low_prec, top_global, ir);
160 default: gmx_incons("Invalid reduced precision file format");
163 if ((EI_DYNAMICS(ir->eI) || EI_ENERGY_MINIMIZATION(ir->eI))
165 && !(EI_DYNAMICS(ir->eI) && ir->nstxout == 0 && ir->nstvout == 0 && ir->nstfout == 0)))
167 const char* filename;
168 filename = ftp2fn(efTRN, nfile, fnm);
169 switch (fn2ftp(filename))
173 /* If there is no uncompressed coordinate output and
174 there is compressed TNG output write forces
175 and/or velocities to the TNG file instead. */
176 if (ir->nstxout != 0 || ir->nstxout_compressed == 0 || !of->tng_low_prec)
178 of->fp_trn = gmx_trr_open(filename, filemode);
182 gmx_tng_open(filename, filemode[0], &of->tng);
183 if (filemode[0] == 'w')
185 gmx_tng_prepare_md_writing(of->tng, top_global, ir);
189 default: gmx_incons("Invalid full precision file format");
192 if (EI_DYNAMICS(ir->eI) || EI_ENERGY_MINIMIZATION(ir->eI))
194 of->fp_ene = open_enx(ftp2fn(efEDR, nfile, fnm), filemode);
196 of->fn_cpt = opt2fn("-cpo", nfile, fnm);
198 if ((ir->efep != efepNO || ir->bSimTemp) && ir->fepvals->nstdhdl > 0
199 && (ir->fepvals->separate_dhdl_file == esepdhdlfileYES) && EI_DYNAMICS(ir->eI))
201 if (restartWithAppending)
203 of->fp_dhdl = gmx_fio_fopen(opt2fn("-dhdl", nfile, fnm), filemode);
207 of->fp_dhdl = open_dhdl(opt2fn("-dhdl", nfile, fnm), ir, oenv);
211 outputProvider->initOutput(fplog, nfile, fnm, restartWithAppending, oenv);
212 of->mdModulesNotifier = &mdModulesNotifier;
214 /* Set up atom counts so they can be passed to actual
215 trajectory-writing routines later. Also, XTC writing needs
216 to know what (and how many) atoms might be in the XTC
217 groups, and how to look up later which ones they are. */
218 of->natoms_global = top_global->natoms;
219 of->groups = &top_global->groups;
220 of->natoms_x_compressed = 0;
221 for (i = 0; (i < top_global->natoms); i++)
223 if (getGroupType(*of->groups, SimulationAtomGroupType::CompressedPositionOutput, i) == 0)
225 of->natoms_x_compressed++;
229 if (ir->nstfout && DOMAINDECOMP(cr))
231 snew(of->f_global, top_global->natoms);
237 please_cite(fplog, "Lundborg2014");
243 ener_file_t mdoutf_get_fp_ene(gmx_mdoutf_t of)
248 FILE* mdoutf_get_fp_dhdl(gmx_mdoutf_t of)
253 gmx_wallcycle_t mdoutf_get_wcycle(gmx_mdoutf_t of)
258 void mdoutf_write_to_trajectory_files(FILE* fplog,
265 t_state* state_local,
266 t_state* state_global,
267 ObservablesHistory* observablesHistory,
268 gmx::ArrayRef<gmx::RVec> f_local)
272 if (DOMAINDECOMP(cr))
274 if (mdof_flags & MDOF_CPT)
276 dd_collect_state(cr->dd, state_local, state_global);
280 if (mdof_flags & (MDOF_X | MDOF_X_COMPRESSED))
282 auto globalXRef = MASTER(cr) ? state_global->x : gmx::ArrayRef<gmx::RVec>();
283 dd_collect_vec(cr->dd, state_local, state_local->x, globalXRef);
285 if (mdof_flags & MDOF_V)
287 auto globalVRef = MASTER(cr) ? state_global->v : gmx::ArrayRef<gmx::RVec>();
288 dd_collect_vec(cr->dd, state_local, state_local->v, globalVRef);
291 f_global = of->f_global;
292 if (mdof_flags & MDOF_F)
294 dd_collect_vec(cr->dd, state_local, f_local,
295 gmx::arrayRefFromArray(reinterpret_cast<gmx::RVec*>(f_global), f_local.size()));
300 /* We have the whole state locally: copy the local state pointer */
301 state_global = state_local;
303 f_global = as_rvec_array(f_local.data());
308 if (mdof_flags & MDOF_CPT)
311 fflush_tng(of->tng_low_prec);
312 /* Write the checkpoint file.
313 * When simulations share the state, an MPI barrier is applied before
314 * renaming old and new checkpoint files to minimize the risk of
315 * checkpoint files getting out of sync.
317 ivec one_ivec = { 1, 1, 1 };
318 write_checkpoint(of->fn_cpt, of->bKeepAndNumCPT, fplog, cr,
319 DOMAINDECOMP(cr) ? cr->dd->numCells : one_ivec,
320 DOMAINDECOMP(cr) ? cr->dd->nnodes : cr->nnodes, of->eIntegrator,
321 of->simulation_part, of->bExpanded, of->elamstats, step, t,
322 state_global, observablesHistory, *(of->mdModulesNotifier),
323 of->simulationsShareState, of->mpiCommMasters);
326 if (mdof_flags & (MDOF_X | MDOF_V | MDOF_F))
328 const rvec* x = (mdof_flags & MDOF_X) ? state_global->x.rvec_array() : nullptr;
329 const rvec* v = (mdof_flags & MDOF_V) ? state_global->v.rvec_array() : nullptr;
330 const rvec* f = (mdof_flags & MDOF_F) ? f_global : nullptr;
334 gmx_trr_write_frame(of->fp_trn, step, t, state_local->lambda[efptFEP],
335 state_local->box, natoms, x, v, f);
336 if (gmx_fio_flush(of->fp_trn) != 0)
338 gmx_file("Cannot write trajectory; maybe you are out of disk space?");
342 /* If a TNG file is open for uncompressed coordinate output also write
343 velocities and forces to it. */
346 gmx_fwrite_tng(of->tng, FALSE, step, t, state_local->lambda[efptFEP],
347 state_local->box, natoms, x, v, f);
349 /* If only a TNG file is open for compressed coordinate output (no uncompressed
350 coordinate output) also write forces and velocities to it. */
351 else if (of->tng_low_prec)
353 gmx_fwrite_tng(of->tng_low_prec, FALSE, step, t, state_local->lambda[efptFEP],
354 state_local->box, natoms, x, v, f);
357 if (mdof_flags & MDOF_X_COMPRESSED)
359 rvec* xxtc = nullptr;
361 if (of->natoms_x_compressed == of->natoms_global)
363 /* We are writing the positions of all of the atoms to
364 the compressed output */
365 xxtc = state_global->x.rvec_array();
369 /* We are writing the positions of only a subset of
370 the atoms to the compressed output, so we have to
371 make a copy of the subset of coordinates. */
374 snew(xxtc, of->natoms_x_compressed);
375 auto x = makeArrayRef(state_global->x);
376 for (i = 0, j = 0; (i < of->natoms_global); i++)
378 if (getGroupType(*of->groups, SimulationAtomGroupType::CompressedPositionOutput, i) == 0)
380 copy_rvec(x[i], xxtc[j++]);
384 if (write_xtc(of->fp_xtc, of->natoms_x_compressed, step, t, state_local->box, xxtc,
385 of->x_compression_precision)
389 "XTC error. This indicates you are out of disk space, or a "
390 "simulation with major instabilities resulting in coordinates "
391 "that are NaN or too large to be represented in the XTC format.\n");
393 gmx_fwrite_tng(of->tng_low_prec, TRUE, step, t, state_local->lambda[efptFEP],
394 state_local->box, of->natoms_x_compressed, xxtc, nullptr, nullptr);
395 if (of->natoms_x_compressed != of->natoms_global)
400 if (mdof_flags & (MDOF_BOX | MDOF_LAMBDA) && !(mdof_flags & (MDOF_X | MDOF_V | MDOF_F)))
406 if (mdof_flags & MDOF_BOX)
408 box = state_local->box;
410 if (mdof_flags & MDOF_LAMBDA)
412 lambda = state_local->lambda[efptFEP];
414 gmx_fwrite_tng(of->tng, FALSE, step, t, lambda, box, natoms, nullptr, nullptr, nullptr);
417 if (mdof_flags & (MDOF_BOX_COMPRESSED | MDOF_LAMBDA_COMPRESSED)
418 && !(mdof_flags & (MDOF_X_COMPRESSED)))
420 if (of->tng_low_prec)
424 if (mdof_flags & MDOF_BOX_COMPRESSED)
426 box = state_local->box;
428 if (mdof_flags & MDOF_LAMBDA_COMPRESSED)
430 lambda = state_local->lambda[efptFEP];
432 gmx_fwrite_tng(of->tng_low_prec, FALSE, step, t, lambda, box, natoms, nullptr,
439 void mdoutf_tng_close(gmx_mdoutf_t of)
441 if (of->tng || of->tng_low_prec)
443 wallcycle_start(of->wcycle, ewcTRAJ);
444 gmx_tng_close(&of->tng);
445 gmx_tng_close(&of->tng_low_prec);
446 wallcycle_stop(of->wcycle, ewcTRAJ);
450 void done_mdoutf(gmx_mdoutf_t of)
452 if (of->fp_ene != nullptr)
454 done_ener_file(of->fp_ene);
458 close_xtc(of->fp_xtc);
462 gmx_trr_close(of->fp_trn);
464 if (of->fp_dhdl != nullptr)
466 gmx_fio_fclose(of->fp_dhdl);
468 of->outputProvider->finishOutput();
469 if (of->f_global != nullptr)
474 gmx_tng_close(&of->tng);
475 gmx_tng_close(&of->tng_low_prec);
480 int mdoutf_get_tng_box_output_interval(gmx_mdoutf_t of)
484 return gmx_tng_get_box_output_interval(of->tng);
489 int mdoutf_get_tng_lambda_output_interval(gmx_mdoutf_t of)
493 return gmx_tng_get_lambda_output_interval(of->tng);
498 int mdoutf_get_tng_compressed_box_output_interval(gmx_mdoutf_t of)
500 if (of->tng_low_prec)
502 return gmx_tng_get_box_output_interval(of->tng_low_prec);
507 int mdoutf_get_tng_compressed_lambda_output_interval(gmx_mdoutf_t of)
509 if (of->tng_low_prec)
511 return gmx_tng_get_lambda_output_interval(of->tng_low_prec);