src/gromacs/fileio/mdoutf.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
   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 #include "mdoutf.h"
  36
  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"
  42 #include "trnio.h"
  43 #include "xtcio.h"
  44 #include "tngio.h"
  45 #include "trajectory_writing.h"
  46 #include "checkpoint.h"
  47 #include "copyrite.h"
  48
  49 #include "gromacs/utility/smalloc.h"
  50
  51 struct gmx_mdoutf {
  52     t_fileio         *fp_trn;
  53     t_fileio         *fp_xtc;
  54     tng_trajectory_t  tng;
  55     tng_trajectory_t  tng_low_prec;
  56     int               x_compression_precision; /* only used by XTC output */
  57     ener_file_t       fp_ene;
  58     const char       *fn_cpt;
  59     gmx_bool          bKeepAndNumCPT;
  60     int               eIntegrator;
  61     gmx_bool          bExpanded;
  62     int               elamstats;
  63     int               simulation_part;
  64     FILE             *fp_dhdl;
  65     FILE             *fp_field;
  66     int               natoms_global;
  67     int               natoms_x_compressed;
  68     gmx_groups_t     *groups; /* for compressed position writing */
  69 };
  70
  71
  72 gmx_mdoutf_t init_mdoutf(FILE *fplog, int nfile, const t_filenm fnm[],
  73                          int mdrun_flags, const t_commrec *cr,
  74                          const t_inputrec *ir, gmx_mtop_t *top_global,
  75                          const output_env_t oenv)
  76 {
  77     gmx_mdoutf_t  of;
  78     char          filemode[3];
  79     gmx_bool      bAppendFiles, bCiteTng = FALSE;
  80     int           i;
  81
  82     snew(of, 1);
  83
  84     of->fp_trn       = NULL;
  85     of->fp_ene       = NULL;
  86     of->fp_xtc       = NULL;
  87     of->tng          = NULL;
  88     of->tng_low_prec = NULL;
  89     of->fp_dhdl      = NULL;
  90     of->fp_field     = NULL;
  91
  92     of->eIntegrator             = ir->eI;
  93     of->bExpanded               = ir->bExpanded;
  94     of->elamstats               = ir->expandedvals->elamstats;
  95     of->simulation_part         = ir->simulation_part;
  96     of->x_compression_precision = ir->x_compression_precision;
  97
  98     if (MASTER(cr))
  99     {
 100         bAppendFiles = (mdrun_flags & MD_APPENDFILES);
 101
 102         of->bKeepAndNumCPT = (mdrun_flags & MD_KEEPANDNUMCPT);
 103
 104         sprintf(filemode, bAppendFiles ? "a+" : "w+");
 105
 106         if ((EI_DYNAMICS(ir->eI) || EI_ENERGY_MINIMIZATION(ir->eI))
 107 #ifndef GMX_FAHCORE
 108             &&
 109             !(EI_DYNAMICS(ir->eI) &&
 110               ir->nstxout == 0 &&
 111               ir->nstvout == 0 &&
 112               ir->nstfout == 0)
 113 #endif
 114             )
 115         {
 116             const char *filename;
 117             filename = ftp2fn(efTRN, nfile, fnm);
 118             switch (fn2ftp(filename))
 119             {
 120                 case efTRR:
 121                 case efTRN:
 122                     of->fp_trn = open_trn(filename, filemode);
 123                     break;
 124                 case efTNG:
 125                     gmx_tng_open(filename, filemode[0], &of->tng);
 126                     if (filemode[0] == 'w')
 127                     {
 128                         gmx_tng_prepare_md_writing(of->tng, top_global, ir);
 129                     }
 130                     bCiteTng = TRUE;
 131                     break;
 132                 default:
 133                     gmx_incons("Invalid full precision file format");
 134             }
 135         }
 136         if (EI_DYNAMICS(ir->eI) &&
 137             ir->nstxout_compressed > 0)
 138         {
 139             const char *filename;
 140             filename = ftp2fn(efCOMPRESSED, nfile, fnm);
 141             switch (fn2ftp(filename))
 142             {
 143                 case efXTC:
 144                     of->fp_xtc                  = open_xtc(filename, filemode);
 145                     break;
 146                 case efTNG:
 147                     gmx_tng_open(filename, filemode[0], &of->tng_low_prec);
 148                     if (filemode[0] == 'w')
 149                     {
 150                         gmx_tng_prepare_low_prec_writing(of->tng_low_prec, top_global, ir);
 151                     }
 152                     bCiteTng = TRUE;
 153                     break;
 154                 default:
 155                     gmx_incons("Invalid reduced precision file format");
 156             }
 157         }
 158         if (EI_DYNAMICS(ir->eI) || EI_ENERGY_MINIMIZATION(ir->eI))
 159         {
 160             of->fp_ene = open_enx(ftp2fn(efEDR, nfile, fnm), filemode);
 161         }
 162         of->fn_cpt = opt2fn("-cpo", nfile, fnm);
 163
 164         if ((ir->efep != efepNO || ir->bSimTemp) && ir->fepvals->nstdhdl > 0 &&
 165             (ir->fepvals->separate_dhdl_file == esepdhdlfileYES ) &&
 166             EI_DYNAMICS(ir->eI))
 167         {
 168             if (bAppendFiles)
 169             {
 170                 of->fp_dhdl = gmx_fio_fopen(opt2fn("-dhdl", nfile, fnm), filemode);
 171             }
 172             else
 173             {
 174                 of->fp_dhdl = open_dhdl(opt2fn("-dhdl", nfile, fnm), ir, oenv);
 175             }
 176         }
 177
 178         if (opt2bSet("-field", nfile, fnm) &&
 179             (ir->ex[XX].n || ir->ex[YY].n || ir->ex[ZZ].n))
 180         {
 181             if (bAppendFiles)
 182             {
 183                 of->fp_dhdl = gmx_fio_fopen(opt2fn("-field", nfile, fnm),
 184                                             filemode);
 185             }
 186             else
 187             {
 188                 of->fp_field = xvgropen(opt2fn("-field", nfile, fnm),
 189                                         "Applied electric field", "Time (ps)",
 190                                         "E (V/nm)", oenv);
 191             }
 192         }
 193
 194         /* Set up atom counts so they can be passed to actual
 195            trajectory-writing routines later. Also, XTC writing needs
 196            to know what (and how many) atoms might be in the XTC
 197            groups, and how to look up later which ones they are. */
 198         of->natoms_global       = top_global->natoms;
 199         of->groups              = &top_global->groups;
 200         of->natoms_x_compressed = 0;
 201         for (i = 0; (i < top_global->natoms); i++)
 202         {
 203             if (ggrpnr(of->groups, egcCompressedX, i) == 0)
 204             {
 205                 of->natoms_x_compressed++;
 206             }
 207         }
 208     }
 209
 210     if (bCiteTng)
 211     {
 212         please_cite(fplog, "Lundborg2014");
 213     }
 214
 215     return of;
 216 }
 217
 218 FILE *mdoutf_get_fp_field(gmx_mdoutf_t of)
 219 {
 220     return of->fp_field;
 221 }
 222
 223 ener_file_t mdoutf_get_fp_ene(gmx_mdoutf_t of)
 224 {
 225     return of->fp_ene;
 226 }
 227
 228 FILE *mdoutf_get_fp_dhdl(gmx_mdoutf_t of)
 229 {
 230     return of->fp_dhdl;
 231 }
 232
 233 void mdoutf_write_to_trajectory_files(FILE *fplog, t_commrec *cr,
 234                                       gmx_mdoutf_t of,
 235                                       int mdof_flags,
 236                                       gmx_mtop_t *top_global,
 237                                       gmx_int64_t step, double t,
 238                                       t_state *state_local, t_state *state_global,
 239                                       rvec *f_local, rvec *f_global)
 240 {
 241     rvec *local_v;
 242     rvec *global_v;
 243
 244     /* MRS -- defining these variables is to manage the difference
 245      * between half step and full step velocities, but there must be a better way . . . */
 246
 247     local_v  = state_local->v;
 248     global_v = state_global->v;
 249
 250     if (DOMAINDECOMP(cr))
 251     {
 252         if (mdof_flags & MDOF_CPT)
 253         {
 254             dd_collect_state(cr->dd, state_local, state_global);
 255         }
 256         else
 257         {
 258             if (mdof_flags & (MDOF_X | MDOF_X_COMPRESSED))
 259             {
 260                 dd_collect_vec(cr->dd, state_local, state_local->x,
 261                                state_global->x);
 262             }
 263             if (mdof_flags & MDOF_V)
 264             {
 265                 dd_collect_vec(cr->dd, state_local, local_v,
 266                                global_v);
 267             }
 268         }
 269         if (mdof_flags & MDOF_F)
 270         {
 271             dd_collect_vec(cr->dd, state_local, f_local, f_global);
 272         }
 273     }
 274     else
 275     {
 276         if (mdof_flags & MDOF_CPT)
 277         {
 278             /* All pointers in state_local are equal to state_global,
 279              * but we need to copy the non-pointer entries.
 280              */
 281             state_global->lambda = state_local->lambda;
 282             state_global->veta   = state_local->veta;
 283             state_global->vol0   = state_local->vol0;
 284             copy_mat(state_local->box, state_global->box);
 285             copy_mat(state_local->boxv, state_global->boxv);
 286             copy_mat(state_local->svir_prev, state_global->svir_prev);
 287             copy_mat(state_local->fvir_prev, state_global->fvir_prev);
 288             copy_mat(state_local->pres_prev, state_global->pres_prev);
 289         }
 290     }
 291
 292     if (MASTER(cr))
 293     {
 294         if (mdof_flags & MDOF_CPT)
 295         {
 296             fflush_tng(of->tng);
 297             fflush_tng(of->tng_low_prec);
 298             write_checkpoint(of->fn_cpt, of->bKeepAndNumCPT,
 299                              fplog, cr, of->eIntegrator, of->simulation_part,
 300                              of->bExpanded, of->elamstats, step, t, state_global);
 301         }
 302
 303         if (mdof_flags & (MDOF_X | MDOF_V | MDOF_F))
 304         {
 305             if (of->fp_trn)
 306             {
 307                 fwrite_trn(of->fp_trn, step, t, state_local->lambda[efptFEP],
 308                            state_local->box, top_global->natoms,
 309                            (mdof_flags & MDOF_X) ? state_global->x : NULL,
 310                            (mdof_flags & MDOF_V) ? global_v : NULL,
 311                            (mdof_flags & MDOF_F) ? f_global : NULL);
 312                 if (gmx_fio_flush(of->fp_trn) != 0)
 313                 {
 314                     gmx_file("Cannot write trajectory; maybe you are out of disk space?");
 315                 }
 316             }
 317
 318             gmx_fwrite_tng(of->tng, FALSE, step, t, state_local->lambda[efptFEP],
 319                            (const rvec *) state_local->box,
 320                            top_global->natoms,
 321                            (mdof_flags & MDOF_X) ? (const rvec *) state_global->x : NULL,
 322                            (mdof_flags & MDOF_V) ? (const rvec *) global_v : NULL,
 323                            (mdof_flags & MDOF_F) ? (const rvec *) f_global : NULL);
 324         }
 325         if (mdof_flags & MDOF_X_COMPRESSED)
 326         {
 327             rvec *xxtc = NULL;
 328
 329             if (of->natoms_x_compressed == of->natoms_global)
 330             {
 331                 /* We are writing the positions of all of the atoms to
 332                    the compressed output */
 333                 xxtc = state_global->x;
 334             }
 335             else
 336             {
 337                 /* We are writing the positions of only a subset of
 338                    the atoms to the compressed output, so we have to
 339                    make a copy of the subset of coordinates. */
 340                 int i, j;
 341
 342                 snew(xxtc, of->natoms_x_compressed);
 343                 for (i = 0, j = 0; (i < of->natoms_x_compressed); i++)
 344                 {
 345                     if (ggrpnr(of->groups, egcCompressedX, i) == 0)
 346                     {
 347                         copy_rvec(state_global->x[i], xxtc[j++]);
 348                     }
 349                 }
 350             }
 351             if (write_xtc(of->fp_xtc, of->natoms_x_compressed, step, t,
 352                           state_local->box, xxtc, of->x_compression_precision) == 0)
 353             {
 354                 gmx_fatal(FARGS, "XTC error - maybe you are out of disk space?");
 355             }
 356             gmx_fwrite_tng(of->tng_low_prec,
 357                            TRUE,
 358                            step,
 359                            t,
 360                            state_local->lambda[efptFEP],
 361                            (const rvec *) state_local->box,
 362                            of->natoms_x_compressed,
 363                            (const rvec *) xxtc,
 364                            NULL,
 365                            NULL);
 366             if (of->natoms_x_compressed != of->natoms_global)
 367             {
 368                 sfree(xxtc);
 369             }
 370         }
 371     }
 372 }
 373
 374 void done_mdoutf(gmx_mdoutf_t of)
 375 {
 376     if (of->fp_ene != NULL)
 377     {
 378         close_enx(of->fp_ene);
 379     }
 380     if (of->fp_xtc)
 381     {
 382         close_xtc(of->fp_xtc);
 383     }
 384     if (of->fp_trn)
 385     {
 386         close_trn(of->fp_trn);
 387     }
 388     if (of->fp_dhdl != NULL)
 389     {
 390         gmx_fio_fclose(of->fp_dhdl);
 391     }
 392     if (of->fp_field != NULL)
 393     {
 394         gmx_fio_fclose(of->fp_field);
 395     }
 396     gmx_tng_close(&of->tng);
 397     gmx_tng_close(&of->tng_low_prec);
 398
 399     sfree(of);
 400 }