src/gromacs/mdlib/mdoutf.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2013,2014,2015,2016, 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 "gmxpre.h"
  36
  37 #include "mdoutf.h"
  38
  39 #include "gromacs/commandline/filenm.h"
  40 #include "gromacs/domdec/domdec.h"
  41 #include "gromacs/domdec/domdec_struct.h"
  42 #include "gromacs/fileio/checkpoint.h"
  43 #include "gromacs/fileio/gmxfio.h"
  44 #include "gromacs/fileio/tngio.h"
  45 #include "gromacs/fileio/trrio.h"
  46 #include "gromacs/fileio/xtcio.h"
  47 #include "gromacs/fileio/xvgr.h"
  48 #include "gromacs/math/vec.h"
  49 #include "gromacs/mdlib/mdrun.h"
  50 #include "gromacs/mdlib/trajectory_writing.h"
  51 #include "gromacs/mdtypes/commrec.h"
  52 #include "gromacs/mdtypes/inputrec.h"
  53 #include "gromacs/mdtypes/md_enums.h"
  54 #include "gromacs/timing/wallcycle.h"
  55 #include "gromacs/utility/fatalerror.h"
  56 #include "gromacs/utility/pleasecite.h"
  57 #include "gromacs/utility/smalloc.h"
  58
  59 struct gmx_mdoutf {
  60     t_fileio         *fp_trn;
  61     t_fileio         *fp_xtc;
  62     tng_trajectory_t  tng;
  63     tng_trajectory_t  tng_low_prec;
  64     int               x_compression_precision; /* only used by XTC output */
  65     ener_file_t       fp_ene;
  66     const char       *fn_cpt;
  67     gmx_bool          bKeepAndNumCPT;
  68     int               eIntegrator;
  69     gmx_bool          bExpanded;
  70     int               elamstats;
  71     int               simulation_part;
  72     FILE             *fp_dhdl;
  73     int               natoms_global;
  74     int               natoms_x_compressed;
  75     gmx_groups_t     *groups; /* for compressed position writing */
  76     gmx_wallcycle_t   wcycle;
  77     rvec             *f_global;
  78 };
  79
  80
  81 gmx_mdoutf_t init_mdoutf(FILE *fplog, int nfile, const t_filenm fnm[],
  82                          int mdrun_flags, const t_commrec *cr,
  83                          const t_inputrec *ir, gmx_mtop_t *top_global,
  84                          const gmx_output_env_t *oenv, gmx_wallcycle_t wcycle)
  85 {
  86     gmx_mdoutf_t   of;
  87     const char    *appendMode = "a+", *writeMode = "w+", *filemode;
  88     gmx_bool       bAppendFiles, bCiteTng = FALSE;
  89     int            i;
  90
  91     snew(of, 1);
  92
  93     of->fp_trn       = NULL;
  94     of->fp_ene       = NULL;
  95     of->fp_xtc       = NULL;
  96     of->tng          = NULL;
  97     of->tng_low_prec = NULL;
  98     of->fp_dhdl      = NULL;
  99
 100     of->eIntegrator             = ir->eI;
 101     of->bExpanded               = ir->bExpanded;
 102     of->elamstats               = ir->expandedvals->elamstats;
 103     of->simulation_part         = ir->simulation_part;
 104     of->x_compression_precision = static_cast<int>(ir->x_compression_precision);
 105     of->wcycle                  = wcycle;
 106     of->f_global                = NULL;
 107
 108     if (MASTER(cr))
 109     {
 110         bAppendFiles = (mdrun_flags & MD_APPENDFILES);
 111
 112         of->bKeepAndNumCPT = (mdrun_flags & MD_KEEPANDNUMCPT);
 113
 114         filemode = bAppendFiles ? appendMode : writeMode;
 115
 116         if (EI_DYNAMICS(ir->eI) &&
 117             ir->nstxout_compressed > 0)
 118         {
 119             const char *filename;
 120             filename = ftp2fn(efCOMPRESSED, nfile, fnm);
 121             switch (fn2ftp(filename))
 122             {
 123                 case efXTC:
 124                     of->fp_xtc                  = open_xtc(filename, filemode);
 125                     break;
 126                 case efTNG:
 127                     gmx_tng_open(filename, filemode[0], &of->tng_low_prec);
 128                     if (filemode[0] == 'w')
 129                     {
 130                         gmx_tng_prepare_low_prec_writing(of->tng_low_prec, top_global, ir);
 131                     }
 132                     bCiteTng = TRUE;
 133                     break;
 134                 default:
 135                     gmx_incons("Invalid reduced precision file format");
 136             }
 137         }
 138         if ((EI_DYNAMICS(ir->eI) || EI_ENERGY_MINIMIZATION(ir->eI))
 139 #ifndef GMX_FAHCORE
 140             &&
 141             !(EI_DYNAMICS(ir->eI) &&
 142               ir->nstxout == 0 &&
 143               ir->nstvout == 0 &&
 144               ir->nstfout == 0)
 145 #endif
 146             )
 147         {
 148             const char *filename;
 149             filename = ftp2fn(efTRN, nfile, fnm);
 150             switch (fn2ftp(filename))
 151             {
 152                 case efTRR:
 153                 case efTRN:
 154                     /* If there is no uncompressed coordinate output and
 155                        there is compressed TNG output write forces
 156                        and/or velocities to the TNG file instead. */
 157                     if (ir->nstxout != 0 || ir->nstxout_compressed == 0 ||
 158                         !of->tng_low_prec)
 159                     {
 160                         of->fp_trn = gmx_trr_open(filename, filemode);
 161                     }
 162                     break;
 163                 case efTNG:
 164                     gmx_tng_open(filename, filemode[0], &of->tng);
 165                     if (filemode[0] == 'w')
 166                     {
 167                         gmx_tng_prepare_md_writing(of->tng, top_global, ir);
 168                     }
 169                     bCiteTng = TRUE;
 170                     break;
 171                 default:
 172                     gmx_incons("Invalid full precision file format");
 173             }
 174         }
 175         if (EI_DYNAMICS(ir->eI) || EI_ENERGY_MINIMIZATION(ir->eI))
 176         {
 177             of->fp_ene = open_enx(ftp2fn(efEDR, nfile, fnm), filemode);
 178         }
 179         of->fn_cpt = opt2fn("-cpo", nfile, fnm);
 180
 181         if ((ir->efep != efepNO || ir->bSimTemp) && ir->fepvals->nstdhdl > 0 &&
 182             (ir->fepvals->separate_dhdl_file == esepdhdlfileYES ) &&
 183             EI_DYNAMICS(ir->eI))
 184         {
 185             if (bAppendFiles)
 186             {
 187                 of->fp_dhdl = gmx_fio_fopen(opt2fn("-dhdl", nfile, fnm), filemode);
 188             }
 189             else
 190             {
 191                 of->fp_dhdl = open_dhdl(opt2fn("-dhdl", nfile, fnm), ir, oenv);
 192             }
 193         }
 194
 195         ir->efield->initOutput(fplog, nfile, fnm, bAppendFiles, oenv);
 196
 197         /* Set up atom counts so they can be passed to actual
 198            trajectory-writing routines later. Also, XTC writing needs
 199            to know what (and how many) atoms might be in the XTC
 200            groups, and how to look up later which ones they are. */
 201         of->natoms_global       = top_global->natoms;
 202         of->groups              = &top_global->groups;
 203         of->natoms_x_compressed = 0;
 204         for (i = 0; (i < top_global->natoms); i++)
 205         {
 206             if (ggrpnr(of->groups, egcCompressedX, i) == 0)
 207             {
 208                 of->natoms_x_compressed++;
 209             }
 210         }
 211
 212         if (ir->nstfout && DOMAINDECOMP(cr))
 213         {
 214             snew(of->f_global, top_global->natoms);
 215         }
 216     }
 217
 218     if (bCiteTng)
 219     {
 220         please_cite(fplog, "Lundborg2014");
 221     }
 222
 223     return of;
 224 }
 225
 226 ener_file_t mdoutf_get_fp_ene(gmx_mdoutf_t of)
 227 {
 228     return of->fp_ene;
 229 }
 230
 231 FILE *mdoutf_get_fp_dhdl(gmx_mdoutf_t of)
 232 {
 233     return of->fp_dhdl;
 234 }
 235
 236 gmx_wallcycle_t mdoutf_get_wcycle(gmx_mdoutf_t of)
 237 {
 238     return of->wcycle;
 239 }
 240
 241 void mdoutf_write_to_trajectory_files(FILE *fplog, t_commrec *cr,
 242                                       gmx_mdoutf_t of,
 243                                       int mdof_flags,
 244                                       gmx_mtop_t *top_global,
 245                                       gmx_int64_t step, double t,
 246                                       t_state *state_local, t_state *state_global,
 247                                       energyhistory_t *energyHistory,
 248                                       PaddedRVecVector *f_local)
 249 {
 250     rvec *f_global;
 251
 252     if (DOMAINDECOMP(cr))
 253     {
 254         if (mdof_flags & MDOF_CPT)
 255         {
 256             dd_collect_state(cr->dd, state_local, state_global);
 257         }
 258         else
 259         {
 260             if (mdof_flags & (MDOF_X | MDOF_X_COMPRESSED))
 261             {
 262                 dd_collect_vec(cr->dd, state_local, &state_local->x,
 263                                &state_global->x);
 264             }
 265             if (mdof_flags & MDOF_V)
 266             {
 267                 dd_collect_vec(cr->dd, state_local, &state_local->v,
 268                                &state_global->v);
 269             }
 270         }
 271         f_global = of->f_global;
 272         if (mdof_flags & MDOF_F)
 273         {
 274             dd_collect_vec(cr->dd, state_local, f_local, f_global);
 275         }
 276     }
 277     else
 278     {
 279         /* We have the whole state locally: copy the local state pointer */
 280         state_global = state_local;
 281
 282         f_global     = as_rvec_array(f_local->data());
 283     }
 284
 285     if (MASTER(cr))
 286     {
 287         if (mdof_flags & MDOF_CPT)
 288         {
 289             fflush_tng(of->tng);
 290             fflush_tng(of->tng_low_prec);
 291             ivec one_ivec = { 1, 1, 1 };
 292             write_checkpoint(of->fn_cpt, of->bKeepAndNumCPT,
 293                              fplog, cr,
 294                              DOMAINDECOMP(cr) ? cr->dd->nc : one_ivec,
 295                              DOMAINDECOMP(cr) ? cr->dd->nnodes : cr->nnodes,
 296                              of->eIntegrator, of->simulation_part,
 297                              of->bExpanded, of->elamstats, step, t,
 298                              state_global, energyHistory);
 299         }
 300
 301         if (mdof_flags & (MDOF_X | MDOF_V | MDOF_F))
 302         {
 303             const rvec *x = (mdof_flags & MDOF_X) ? as_rvec_array(state_global->x.data()) : NULL;
 304             const rvec *v = (mdof_flags & MDOF_V) ? as_rvec_array(state_global->v.data()) : NULL;
 305             const rvec *f = (mdof_flags & MDOF_F) ? f_global : NULL;
 306
 307             if (of->fp_trn)
 308             {
 309                 gmx_trr_write_frame(of->fp_trn, step, t, state_local->lambda[efptFEP],
 310                                     state_local->box, top_global->natoms,
 311                                     x, v, f);
 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             /* If a TNG file is open for uncompressed coordinate output also write
 319                velocities and forces to it. */
 320             else if (of->tng)
 321             {
 322                 gmx_fwrite_tng(of->tng, FALSE, step, t, state_local->lambda[efptFEP],
 323                                state_local->box,
 324                                top_global->natoms,
 325                                x, v, f);
 326             }
 327             /* If only a TNG file is open for compressed coordinate output (no uncompressed
 328                coordinate output) also write forces and velocities to it. */
 329             else if (of->tng_low_prec)
 330             {
 331                 gmx_fwrite_tng(of->tng_low_prec, FALSE, step, t, state_local->lambda[efptFEP],
 332                                state_local->box,
 333                                top_global->natoms,
 334                                x, v, f);
 335             }
 336         }
 337         if (mdof_flags & MDOF_X_COMPRESSED)
 338         {
 339             rvec *xxtc = NULL;
 340
 341             if (of->natoms_x_compressed == of->natoms_global)
 342             {
 343                 /* We are writing the positions of all of the atoms to
 344                    the compressed output */
 345                 xxtc = as_rvec_array(state_global->x.data());
 346             }
 347             else
 348             {
 349                 /* We are writing the positions of only a subset of
 350                    the atoms to the compressed output, so we have to
 351                    make a copy of the subset of coordinates. */
 352                 int i, j;
 353
 354                 snew(xxtc, of->natoms_x_compressed);
 355                 for (i = 0, j = 0; (i < of->natoms_global); i++)
 356                 {
 357                     if (ggrpnr(of->groups, egcCompressedX, i) == 0)
 358                     {
 359                         copy_rvec(state_global->x[i], xxtc[j++]);
 360                     }
 361                 }
 362             }
 363             if (write_xtc(of->fp_xtc, of->natoms_x_compressed, step, t,
 364                           state_local->box, xxtc, of->x_compression_precision) == 0)
 365             {
 366                 gmx_fatal(FARGS, "XTC error - maybe you are out of disk space?");
 367             }
 368             gmx_fwrite_tng(of->tng_low_prec,
 369                            TRUE,
 370                            step,
 371                            t,
 372                            state_local->lambda[efptFEP],
 373                            state_local->box,
 374                            of->natoms_x_compressed,
 375                            xxtc,
 376                            NULL,
 377                            NULL);
 378             if (of->natoms_x_compressed != of->natoms_global)
 379             {
 380                 sfree(xxtc);
 381             }
 382         }
 383     }
 384 }
 385
 386 void mdoutf_tng_close(gmx_mdoutf_t of)
 387 {
 388     if (of->tng || of->tng_low_prec)
 389     {
 390         wallcycle_start(of->wcycle, ewcTRAJ);
 391         gmx_tng_close(&of->tng);
 392         gmx_tng_close(&of->tng_low_prec);
 393         wallcycle_stop(of->wcycle, ewcTRAJ);
 394     }
 395 }
 396
 397 void done_mdoutf(gmx_mdoutf_t of, const t_inputrec *ir)
 398 {
 399     if (of->fp_ene != NULL)
 400     {
 401         close_enx(of->fp_ene);
 402     }
 403     if (of->fp_xtc)
 404     {
 405         close_xtc(of->fp_xtc);
 406     }
 407     if (of->fp_trn)
 408     {
 409         gmx_trr_close(of->fp_trn);
 410     }
 411     if (of->fp_dhdl != NULL)
 412     {
 413         gmx_fio_fclose(of->fp_dhdl);
 414     }
 415     ir->efield->finishOutput();
 416     if (of->f_global != NULL)
 417     {
 418         sfree(of->f_global);
 419     }
 420
 421     gmx_tng_close(&of->tng);
 422     gmx_tng_close(&of->tng_low_prec);
 423
 424     sfree(of);
 425 }