2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2013,2014, 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 #include "gromacs/legacyheaders/xvgr.h"
38 #include "gromacs/legacyheaders/mdrun.h"
39 #include "gromacs/legacyheaders/types/commrec.h"
40 #include "gromacs/legacyheaders/mvdata.h"
41 #include "gromacs/legacyheaders/domdec.h"
45 #include "trajectory_writing.h"
46 #include "checkpoint.h"
49 #include "gromacs/legacyheaders/gmx_fatal.h"
50 #include "gromacs/utility/smalloc.h"
51 #include "gromacs/timing/wallcycle.h"
57 tng_trajectory_t tng_low_prec;
58 int x_compression_precision; /* only used by XTC output */
61 gmx_bool bKeepAndNumCPT;
69 int natoms_x_compressed;
70 gmx_groups_t *groups; /* for compressed position writing */
71 gmx_wallcycle_t wcycle;
75 gmx_mdoutf_t init_mdoutf(FILE *fplog, int nfile, const t_filenm fnm[],
76 int mdrun_flags, const t_commrec *cr,
77 const t_inputrec *ir, gmx_mtop_t *top_global,
78 const output_env_t oenv, gmx_wallcycle_t wcycle)
82 gmx_bool bAppendFiles, bCiteTng = FALSE;
91 of->tng_low_prec = NULL;
95 of->eIntegrator = ir->eI;
96 of->bExpanded = ir->bExpanded;
97 of->elamstats = ir->expandedvals->elamstats;
98 of->simulation_part = ir->simulation_part;
99 of->x_compression_precision = ir->x_compression_precision;
104 bAppendFiles = (mdrun_flags & MD_APPENDFILES);
106 of->bKeepAndNumCPT = (mdrun_flags & MD_KEEPANDNUMCPT);
108 sprintf(filemode, bAppendFiles ? "a+" : "w+");
110 if ((EI_DYNAMICS(ir->eI) || EI_ENERGY_MINIMIZATION(ir->eI))
113 !(EI_DYNAMICS(ir->eI) &&
120 const char *filename;
121 filename = ftp2fn(efTRN, nfile, fnm);
122 switch (fn2ftp(filename))
126 of->fp_trn = open_trn(filename, filemode);
129 gmx_tng_open(filename, filemode[0], &of->tng);
130 if (filemode[0] == 'w')
132 gmx_tng_prepare_md_writing(of->tng, top_global, ir);
137 gmx_incons("Invalid full precision file format");
140 if (EI_DYNAMICS(ir->eI) &&
141 ir->nstxout_compressed > 0)
143 const char *filename;
144 filename = ftp2fn(efCOMPRESSED, nfile, fnm);
145 switch (fn2ftp(filename))
148 of->fp_xtc = open_xtc(filename, filemode);
151 gmx_tng_open(filename, filemode[0], &of->tng_low_prec);
152 if (filemode[0] == 'w')
154 gmx_tng_prepare_low_prec_writing(of->tng_low_prec, top_global, ir);
159 gmx_incons("Invalid reduced precision file format");
162 if (EI_DYNAMICS(ir->eI) || EI_ENERGY_MINIMIZATION(ir->eI))
164 of->fp_ene = open_enx(ftp2fn(efEDR, nfile, fnm), filemode);
166 of->fn_cpt = opt2fn("-cpo", nfile, fnm);
168 if ((ir->efep != efepNO || ir->bSimTemp) && ir->fepvals->nstdhdl > 0 &&
169 (ir->fepvals->separate_dhdl_file == esepdhdlfileYES ) &&
174 of->fp_dhdl = gmx_fio_fopen(opt2fn("-dhdl", nfile, fnm), filemode);
178 of->fp_dhdl = open_dhdl(opt2fn("-dhdl", nfile, fnm), ir, oenv);
182 if (opt2bSet("-field", nfile, fnm) &&
183 (ir->ex[XX].n || ir->ex[YY].n || ir->ex[ZZ].n))
187 of->fp_dhdl = gmx_fio_fopen(opt2fn("-field", nfile, fnm),
192 of->fp_field = xvgropen(opt2fn("-field", nfile, fnm),
193 "Applied electric field", "Time (ps)",
198 /* Set up atom counts so they can be passed to actual
199 trajectory-writing routines later. Also, XTC writing needs
200 to know what (and how many) atoms might be in the XTC
201 groups, and how to look up later which ones they are. */
202 of->natoms_global = top_global->natoms;
203 of->groups = &top_global->groups;
204 of->natoms_x_compressed = 0;
205 for (i = 0; (i < top_global->natoms); i++)
207 if (ggrpnr(of->groups, egcCompressedX, i) == 0)
209 of->natoms_x_compressed++;
216 please_cite(fplog, "Lundborg2014");
222 FILE *mdoutf_get_fp_field(gmx_mdoutf_t of)
227 ener_file_t mdoutf_get_fp_ene(gmx_mdoutf_t of)
232 FILE *mdoutf_get_fp_dhdl(gmx_mdoutf_t of)
237 gmx_wallcycle_t mdoutf_get_wcycle(gmx_mdoutf_t of)
242 void mdoutf_write_to_trajectory_files(FILE *fplog, t_commrec *cr,
245 gmx_mtop_t *top_global,
246 gmx_int64_t step, double t,
247 t_state *state_local, t_state *state_global,
248 rvec *f_local, rvec *f_global)
253 /* MRS -- defining these variables is to manage the difference
254 * between half step and full step velocities, but there must be a better way . . . */
256 local_v = state_local->v;
257 global_v = state_global->v;
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 dd_collect_vec(cr->dd, state_local, state_local->x,
272 if (mdof_flags & MDOF_V)
274 dd_collect_vec(cr->dd, state_local, local_v,
278 if (mdof_flags & MDOF_F)
280 dd_collect_vec(cr->dd, state_local, f_local, f_global);
285 if (mdof_flags & MDOF_CPT)
287 /* All pointers in state_local are equal to state_global,
288 * but we need to copy the non-pointer entries.
290 state_global->lambda = state_local->lambda;
291 state_global->veta = state_local->veta;
292 state_global->vol0 = state_local->vol0;
293 copy_mat(state_local->box, state_global->box);
294 copy_mat(state_local->boxv, state_global->boxv);
295 copy_mat(state_local->svir_prev, state_global->svir_prev);
296 copy_mat(state_local->fvir_prev, state_global->fvir_prev);
297 copy_mat(state_local->pres_prev, state_global->pres_prev);
303 if (mdof_flags & MDOF_CPT)
306 fflush_tng(of->tng_low_prec);
307 write_checkpoint(of->fn_cpt, of->bKeepAndNumCPT,
308 fplog, cr, of->eIntegrator, of->simulation_part,
309 of->bExpanded, of->elamstats, step, t, state_global);
312 if (mdof_flags & (MDOF_X | MDOF_V | MDOF_F))
316 fwrite_trn(of->fp_trn, step, t, state_local->lambda[efptFEP],
317 state_local->box, top_global->natoms,
318 (mdof_flags & MDOF_X) ? state_global->x : NULL,
319 (mdof_flags & MDOF_V) ? global_v : NULL,
320 (mdof_flags & MDOF_F) ? f_global : NULL);
321 if (gmx_fio_flush(of->fp_trn) != 0)
323 gmx_file("Cannot write trajectory; maybe you are out of disk space?");
327 gmx_fwrite_tng(of->tng, FALSE, step, t, state_local->lambda[efptFEP],
328 (const rvec *) state_local->box,
330 (mdof_flags & MDOF_X) ? (const rvec *) state_global->x : NULL,
331 (mdof_flags & MDOF_V) ? (const rvec *) global_v : NULL,
332 (mdof_flags & MDOF_F) ? (const rvec *) f_global : NULL);
334 if (mdof_flags & MDOF_X_COMPRESSED)
338 if (of->natoms_x_compressed == of->natoms_global)
340 /* We are writing the positions of all of the atoms to
341 the compressed output */
342 xxtc = state_global->x;
346 /* We are writing the positions of only a subset of
347 the atoms to the compressed output, so we have to
348 make a copy of the subset of coordinates. */
351 snew(xxtc, of->natoms_x_compressed);
352 for (i = 0, j = 0; (i < of->natoms_global); i++)
354 if (ggrpnr(of->groups, egcCompressedX, i) == 0)
356 copy_rvec(state_global->x[i], xxtc[j++]);
360 if (write_xtc(of->fp_xtc, of->natoms_x_compressed, step, t,
361 state_local->box, xxtc, of->x_compression_precision) == 0)
363 gmx_fatal(FARGS, "XTC error - maybe you are out of disk space?");
365 gmx_fwrite_tng(of->tng_low_prec,
369 state_local->lambda[efptFEP],
370 (const rvec *) state_local->box,
371 of->natoms_x_compressed,
375 if (of->natoms_x_compressed != of->natoms_global)
383 void mdoutf_tng_close(gmx_mdoutf_t of)
385 if (of->tng || of->tng_low_prec)
387 wallcycle_start(of->wcycle, ewcTRAJ);
388 gmx_tng_close(&of->tng);
389 gmx_tng_close(&of->tng_low_prec);
390 wallcycle_stop(of->wcycle, ewcTRAJ);
394 void done_mdoutf(gmx_mdoutf_t of)
396 if (of->fp_ene != NULL)
398 close_enx(of->fp_ene);
402 close_xtc(of->fp_xtc);
406 close_trn(of->fp_trn);
408 if (of->fp_dhdl != NULL)
410 gmx_fio_fclose(of->fp_dhdl);
412 if (of->fp_field != NULL)
414 gmx_fio_fclose(of->fp_field);
417 gmx_tng_close(&of->tng);
418 gmx_tng_close(&of->tng_low_prec);