2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2013,2014,2015,2016,2017 by 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/mdtypes/commrec.h"
53 #include "gromacs/mdtypes/imdoutputprovider.h"
54 #include "gromacs/mdtypes/inputrec.h"
55 #include "gromacs/mdtypes/md_enums.h"
56 #include "gromacs/mdtypes/mdrunoptions.h"
57 #include "gromacs/mdtypes/state.h"
58 #include "gromacs/timing/wallcycle.h"
59 #include "gromacs/topology/topology.h"
60 #include "gromacs/utility/fatalerror.h"
61 #include "gromacs/utility/pleasecite.h"
62 #include "gromacs/utility/smalloc.h"
68 gmx_tng_trajectory_t tng;
69 gmx_tng_trajectory_t tng_low_prec;
70 int x_compression_precision; /* only used by XTC output */
73 gmx_bool bKeepAndNumCPT;
80 int natoms_x_compressed;
81 SimulationGroups* groups; /* for compressed position writing */
82 gmx_wallcycle_t wcycle;
84 gmx::IMDOutputProvider* outputProvider;
85 const gmx::MdModulesNotifier* mdModulesNotifier;
89 gmx_mdoutf_t init_mdoutf(FILE* fplog,
92 const gmx::MdrunOptions& mdrunOptions,
94 gmx::IMDOutputProvider* outputProvider,
95 const gmx::MdModulesNotifier& mdModulesNotifier,
97 gmx_mtop_t* top_global,
98 const gmx_output_env_t* oenv,
99 gmx_wallcycle_t wcycle,
100 const gmx::StartingBehavior startingBehavior)
103 const char * appendMode = "a+", *writeMode = "w+", *filemode;
104 gmx_bool bCiteTng = FALSE;
106 bool restartWithAppending = (startingBehavior == gmx::StartingBehavior::RestartWithAppending);
110 of->fp_trn = nullptr;
111 of->fp_ene = nullptr;
112 of->fp_xtc = nullptr;
114 of->tng_low_prec = nullptr;
115 of->fp_dhdl = nullptr;
117 of->eIntegrator = ir->eI;
118 of->bExpanded = ir->bExpanded;
119 of->elamstats = ir->expandedvals->elamstats;
120 of->simulation_part = ir->simulation_part;
121 of->x_compression_precision = static_cast<int>(ir->x_compression_precision);
123 of->f_global = nullptr;
124 of->outputProvider = outputProvider;
128 of->bKeepAndNumCPT = mdrunOptions.checkpointOptions.keepAndNumberCheckpointFiles;
130 filemode = restartWithAppending ? appendMode : writeMode;
132 if (EI_DYNAMICS(ir->eI) && ir->nstxout_compressed > 0)
134 const char* filename;
135 filename = ftp2fn(efCOMPRESSED, nfile, fnm);
136 switch (fn2ftp(filename))
138 case efXTC: of->fp_xtc = open_xtc(filename, filemode); break;
140 gmx_tng_open(filename, filemode[0], &of->tng_low_prec);
141 if (filemode[0] == 'w')
143 gmx_tng_prepare_low_prec_writing(of->tng_low_prec, top_global, ir);
147 default: gmx_incons("Invalid reduced precision file format");
150 if ((EI_DYNAMICS(ir->eI) || EI_ENERGY_MINIMIZATION(ir->eI))
152 && !(EI_DYNAMICS(ir->eI) && ir->nstxout == 0 && ir->nstvout == 0 && ir->nstfout == 0)))
154 const char* filename;
155 filename = ftp2fn(efTRN, nfile, fnm);
156 switch (fn2ftp(filename))
160 /* If there is no uncompressed coordinate output and
161 there is compressed TNG output write forces
162 and/or velocities to the TNG file instead. */
163 if (ir->nstxout != 0 || ir->nstxout_compressed == 0 || !of->tng_low_prec)
165 of->fp_trn = gmx_trr_open(filename, filemode);
169 gmx_tng_open(filename, filemode[0], &of->tng);
170 if (filemode[0] == 'w')
172 gmx_tng_prepare_md_writing(of->tng, top_global, ir);
176 default: gmx_incons("Invalid full precision file format");
179 if (EI_DYNAMICS(ir->eI) || EI_ENERGY_MINIMIZATION(ir->eI))
181 of->fp_ene = open_enx(ftp2fn(efEDR, nfile, fnm), filemode);
183 of->fn_cpt = opt2fn("-cpo", nfile, fnm);
185 if ((ir->efep != efepNO || ir->bSimTemp) && ir->fepvals->nstdhdl > 0
186 && (ir->fepvals->separate_dhdl_file == esepdhdlfileYES) && EI_DYNAMICS(ir->eI))
188 if (restartWithAppending)
190 of->fp_dhdl = gmx_fio_fopen(opt2fn("-dhdl", nfile, fnm), filemode);
194 of->fp_dhdl = open_dhdl(opt2fn("-dhdl", nfile, fnm), ir, oenv);
198 outputProvider->initOutput(fplog, nfile, fnm, restartWithAppending, oenv);
199 of->mdModulesNotifier = &mdModulesNotifier;
201 /* Set up atom counts so they can be passed to actual
202 trajectory-writing routines later. Also, XTC writing needs
203 to know what (and how many) atoms might be in the XTC
204 groups, and how to look up later which ones they are. */
205 of->natoms_global = top_global->natoms;
206 of->groups = &top_global->groups;
207 of->natoms_x_compressed = 0;
208 for (i = 0; (i < top_global->natoms); i++)
210 if (getGroupType(*of->groups, SimulationAtomGroupType::CompressedPositionOutput, i) == 0)
212 of->natoms_x_compressed++;
216 if (ir->nstfout && DOMAINDECOMP(cr))
218 snew(of->f_global, top_global->natoms);
224 please_cite(fplog, "Lundborg2014");
230 ener_file_t mdoutf_get_fp_ene(gmx_mdoutf_t of)
235 FILE* mdoutf_get_fp_dhdl(gmx_mdoutf_t of)
240 gmx_wallcycle_t mdoutf_get_wcycle(gmx_mdoutf_t of)
245 void mdoutf_write_to_trajectory_files(FILE* fplog,
252 t_state* state_local,
253 t_state* state_global,
254 ObservablesHistory* observablesHistory,
255 gmx::ArrayRef<gmx::RVec> f_local)
259 if (DOMAINDECOMP(cr))
261 if (mdof_flags & MDOF_CPT)
263 dd_collect_state(cr->dd, state_local, state_global);
267 if (mdof_flags & (MDOF_X | MDOF_X_COMPRESSED))
269 auto globalXRef = MASTER(cr) ? state_global->x : gmx::ArrayRef<gmx::RVec>();
270 dd_collect_vec(cr->dd, state_local, state_local->x, globalXRef);
272 if (mdof_flags & MDOF_V)
274 auto globalVRef = MASTER(cr) ? state_global->v : gmx::ArrayRef<gmx::RVec>();
275 dd_collect_vec(cr->dd, state_local, state_local->v, globalVRef);
278 f_global = of->f_global;
279 if (mdof_flags & MDOF_F)
281 dd_collect_vec(cr->dd, state_local, f_local,
282 gmx::arrayRefFromArray(reinterpret_cast<gmx::RVec*>(f_global), f_local.size()));
287 /* We have the whole state locally: copy the local state pointer */
288 state_global = state_local;
290 f_global = as_rvec_array(f_local.data());
295 if (mdof_flags & MDOF_CPT)
298 fflush_tng(of->tng_low_prec);
299 ivec one_ivec = { 1, 1, 1 };
300 write_checkpoint(of->fn_cpt, of->bKeepAndNumCPT, fplog, cr,
301 DOMAINDECOMP(cr) ? cr->dd->numCells : one_ivec,
302 DOMAINDECOMP(cr) ? cr->dd->nnodes : cr->nnodes, of->eIntegrator,
303 of->simulation_part, of->bExpanded, of->elamstats, step, t,
304 state_global, observablesHistory, *(of->mdModulesNotifier));
307 if (mdof_flags & (MDOF_X | MDOF_V | MDOF_F))
309 const rvec* x = (mdof_flags & MDOF_X) ? state_global->x.rvec_array() : nullptr;
310 const rvec* v = (mdof_flags & MDOF_V) ? state_global->v.rvec_array() : nullptr;
311 const rvec* f = (mdof_flags & MDOF_F) ? f_global : nullptr;
315 gmx_trr_write_frame(of->fp_trn, step, t, state_local->lambda[efptFEP],
316 state_local->box, natoms, x, v, f);
317 if (gmx_fio_flush(of->fp_trn) != 0)
319 gmx_file("Cannot write trajectory; maybe you are out of disk space?");
323 /* If a TNG file is open for uncompressed coordinate output also write
324 velocities and forces to it. */
327 gmx_fwrite_tng(of->tng, FALSE, step, t, state_local->lambda[efptFEP],
328 state_local->box, natoms, x, v, f);
330 /* If only a TNG file is open for compressed coordinate output (no uncompressed
331 coordinate output) also write forces and velocities to it. */
332 else if (of->tng_low_prec)
334 gmx_fwrite_tng(of->tng_low_prec, FALSE, step, t, state_local->lambda[efptFEP],
335 state_local->box, natoms, x, v, f);
338 if (mdof_flags & MDOF_X_COMPRESSED)
340 rvec* xxtc = nullptr;
342 if (of->natoms_x_compressed == of->natoms_global)
344 /* We are writing the positions of all of the atoms to
345 the compressed output */
346 xxtc = state_global->x.rvec_array();
350 /* We are writing the positions of only a subset of
351 the atoms to the compressed output, so we have to
352 make a copy of the subset of coordinates. */
355 snew(xxtc, of->natoms_x_compressed);
356 auto x = makeArrayRef(state_global->x);
357 for (i = 0, j = 0; (i < of->natoms_global); i++)
359 if (getGroupType(*of->groups, SimulationAtomGroupType::CompressedPositionOutput, i) == 0)
361 copy_rvec(x[i], xxtc[j++]);
365 if (write_xtc(of->fp_xtc, of->natoms_x_compressed, step, t, state_local->box, xxtc,
366 of->x_compression_precision)
370 "XTC error. This indicates you are out of disk space, or a "
371 "simulation with major instabilities resulting in coordinates "
372 "that are NaN or too large to be represented in the XTC format.\n");
374 gmx_fwrite_tng(of->tng_low_prec, TRUE, step, t, state_local->lambda[efptFEP],
375 state_local->box, of->natoms_x_compressed, xxtc, nullptr, nullptr);
376 if (of->natoms_x_compressed != of->natoms_global)
381 if (mdof_flags & (MDOF_BOX | MDOF_LAMBDA) && !(mdof_flags & (MDOF_X | MDOF_V | MDOF_F)))
387 if (mdof_flags & MDOF_BOX)
389 box = state_local->box;
391 if (mdof_flags & MDOF_LAMBDA)
393 lambda = state_local->lambda[efptFEP];
395 gmx_fwrite_tng(of->tng, FALSE, step, t, lambda, box, natoms, nullptr, nullptr, nullptr);
398 if (mdof_flags & (MDOF_BOX_COMPRESSED | MDOF_LAMBDA_COMPRESSED)
399 && !(mdof_flags & (MDOF_X_COMPRESSED)))
401 if (of->tng_low_prec)
405 if (mdof_flags & MDOF_BOX_COMPRESSED)
407 box = state_local->box;
409 if (mdof_flags & MDOF_LAMBDA_COMPRESSED)
411 lambda = state_local->lambda[efptFEP];
413 gmx_fwrite_tng(of->tng_low_prec, FALSE, step, t, lambda, box, natoms, nullptr,
420 void mdoutf_tng_close(gmx_mdoutf_t of)
422 if (of->tng || of->tng_low_prec)
424 wallcycle_start(of->wcycle, ewcTRAJ);
425 gmx_tng_close(&of->tng);
426 gmx_tng_close(&of->tng_low_prec);
427 wallcycle_stop(of->wcycle, ewcTRAJ);
431 void done_mdoutf(gmx_mdoutf_t of)
433 if (of->fp_ene != nullptr)
435 done_ener_file(of->fp_ene);
439 close_xtc(of->fp_xtc);
443 gmx_trr_close(of->fp_trn);
445 if (of->fp_dhdl != nullptr)
447 gmx_fio_fclose(of->fp_dhdl);
449 of->outputProvider->finishOutput();
450 if (of->f_global != nullptr)
455 gmx_tng_close(&of->tng);
456 gmx_tng_close(&of->tng_low_prec);
461 int mdoutf_get_tng_box_output_interval(gmx_mdoutf_t of)
465 return gmx_tng_get_box_output_interval(of->tng);
470 int mdoutf_get_tng_lambda_output_interval(gmx_mdoutf_t of)
474 return gmx_tng_get_lambda_output_interval(of->tng);
479 int mdoutf_get_tng_compressed_box_output_interval(gmx_mdoutf_t of)
481 if (of->tng_low_prec)
483 return gmx_tng_get_box_output_interval(of->tng_low_prec);
488 int mdoutf_get_tng_compressed_lambda_output_interval(gmx_mdoutf_t of)
490 if (of->tng_low_prec)
492 return gmx_tng_get_lambda_output_interval(of->tng_low_prec);