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,2017,2018, 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/collect.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/imdoutputprovider.h"
  53 #include "gromacs/mdtypes/inputrec.h"
  54 #include "gromacs/mdtypes/md_enums.h"
  55 #include "gromacs/mdtypes/state.h"
  56 #include "gromacs/timing/wallcycle.h"
  57 #include "gromacs/topology/topology.h"
  58 #include "gromacs/utility/fatalerror.h"
  59 #include "gromacs/utility/pleasecite.h"
  60 #include "gromacs/utility/smalloc.h"
  61
  62 struct gmx_mdoutf {
  63     t_fileio               *fp_trn;
  64     t_fileio               *fp_xtc;
  65     gmx_tng_trajectory_t    tng;
  66     gmx_tng_trajectory_t    tng_low_prec;
  67     int                     x_compression_precision; /* only used by XTC output */
  68     ener_file_t             fp_ene;
  69     const char             *fn_cpt;
  70     gmx_bool                bKeepAndNumCPT;
  71     int                     eIntegrator;
  72     gmx_bool                bExpanded;
  73     int                     elamstats;
  74     int                     simulation_part;
  75     FILE                   *fp_dhdl;
  76     int                     natoms_global;
  77     int                     natoms_x_compressed;
  78     gmx_groups_t           *groups; /* for compressed position writing */
  79     gmx_wallcycle_t         wcycle;
  80     rvec                   *f_global;
  81     gmx::IMDOutputProvider *outputProvider;
  82 };
  83
  84
  85 gmx_mdoutf_t init_mdoutf(FILE *fplog, int nfile, const t_filenm fnm[],
  86                          const MdrunOptions &mdrunOptions,
  87                          const t_commrec *cr,
  88                          gmx::IMDOutputProvider *outputProvider,
  89                          const t_inputrec *ir, gmx_mtop_t *top_global,
  90                          const gmx_output_env_t *oenv, gmx_wallcycle_t wcycle)
  91 {
  92     gmx_mdoutf_t   of;
  93     const char    *appendMode = "a+", *writeMode = "w+", *filemode;
  94     gmx_bool       bAppendFiles, bCiteTng = FALSE;
  95     int            i;
  96
  97     snew(of, 1);
  98
  99     of->fp_trn       = nullptr;
 100     of->fp_ene       = nullptr;
 101     of->fp_xtc       = nullptr;
 102     of->tng          = nullptr;
 103     of->tng_low_prec = nullptr;
 104     of->fp_dhdl      = nullptr;
 105
 106     of->eIntegrator             = ir->eI;
 107     of->bExpanded               = ir->bExpanded;
 108     of->elamstats               = ir->expandedvals->elamstats;
 109     of->simulation_part         = ir->simulation_part;
 110     of->x_compression_precision = static_cast<int>(ir->x_compression_precision);
 111     of->wcycle                  = wcycle;
 112     of->f_global                = nullptr;
 113     of->outputProvider          = outputProvider;
 114
 115     if (MASTER(cr))
 116     {
 117         bAppendFiles = mdrunOptions.continuationOptions.appendFiles;
 118
 119         of->bKeepAndNumCPT = mdrunOptions.checkpointOptions.keepAndNumberCheckpointFiles;
 120
 121         filemode = bAppendFiles ? appendMode : writeMode;
 122
 123         if (EI_DYNAMICS(ir->eI) &&
 124             ir->nstxout_compressed > 0)
 125         {
 126             const char *filename;
 127             filename = ftp2fn(efCOMPRESSED, nfile, fnm);
 128             switch (fn2ftp(filename))
 129             {
 130                 case efXTC:
 131                     of->fp_xtc                  = open_xtc(filename, filemode);
 132                     break;
 133                 case efTNG:
 134                     gmx_tng_open(filename, filemode[0], &of->tng_low_prec);
 135                     if (filemode[0] == 'w')
 136                     {
 137                         gmx_tng_prepare_low_prec_writing(of->tng_low_prec, top_global, ir);
 138                     }
 139                     bCiteTng = TRUE;
 140                     break;
 141                 default:
 142                     gmx_incons("Invalid reduced precision file format");
 143             }
 144         }
 145         if ((EI_DYNAMICS(ir->eI) || EI_ENERGY_MINIMIZATION(ir->eI))
 146 #ifndef GMX_FAHCORE
 147             &&
 148             !(EI_DYNAMICS(ir->eI) &&
 149               ir->nstxout == 0 &&
 150               ir->nstvout == 0 &&
 151               ir->nstfout == 0)
 152 #endif
 153             )
 154         {
 155             const char *filename;
 156             filename = ftp2fn(efTRN, nfile, fnm);
 157             switch (fn2ftp(filename))
 158             {
 159                 case efTRR:
 160                 case efTRN:
 161                     /* If there is no uncompressed coordinate output and
 162                        there is compressed TNG output write forces
 163                        and/or velocities to the TNG file instead. */
 164                     if (ir->nstxout != 0 || ir->nstxout_compressed == 0 ||
 165                         !of->tng_low_prec)
 166                     {
 167                         of->fp_trn = gmx_trr_open(filename, filemode);
 168                     }
 169                     break;
 170                 case efTNG:
 171                     gmx_tng_open(filename, filemode[0], &of->tng);
 172                     if (filemode[0] == 'w')
 173                     {
 174                         gmx_tng_prepare_md_writing(of->tng, top_global, ir);
 175                     }
 176                     bCiteTng = TRUE;
 177                     break;
 178                 default:
 179                     gmx_incons("Invalid full precision file format");
 180             }
 181         }
 182         if (EI_DYNAMICS(ir->eI) || EI_ENERGY_MINIMIZATION(ir->eI))
 183         {
 184             of->fp_ene = open_enx(ftp2fn(efEDR, nfile, fnm), filemode);
 185         }
 186         of->fn_cpt = opt2fn("-cpo", nfile, fnm);
 187
 188         if ((ir->efep != efepNO || ir->bSimTemp) && ir->fepvals->nstdhdl > 0 &&
 189             (ir->fepvals->separate_dhdl_file == esepdhdlfileYES ) &&
 190             EI_DYNAMICS(ir->eI))
 191         {
 192             if (bAppendFiles)
 193             {
 194                 of->fp_dhdl = gmx_fio_fopen(opt2fn("-dhdl", nfile, fnm), filemode);
 195             }
 196             else
 197             {
 198                 of->fp_dhdl = open_dhdl(opt2fn("-dhdl", nfile, fnm), ir, oenv);
 199             }
 200         }
 201
 202         outputProvider->initOutput(fplog, nfile, fnm, bAppendFiles, oenv);
 203
 204         /* Set up atom counts so they can be passed to actual
 205            trajectory-writing routines later. Also, XTC writing needs
 206            to know what (and how many) atoms might be in the XTC
 207            groups, and how to look up later which ones they are. */
 208         of->natoms_global       = top_global->natoms;
 209         of->groups              = &top_global->groups;
 210         of->natoms_x_compressed = 0;
 211         for (i = 0; (i < top_global->natoms); i++)
 212         {
 213             if (ggrpnr(of->groups, egcCompressedX, i) == 0)
 214             {
 215                 of->natoms_x_compressed++;
 216             }
 217         }
 218
 219         if (ir->nstfout && DOMAINDECOMP(cr))
 220         {
 221             snew(of->f_global, top_global->natoms);
 222         }
 223     }
 224
 225     if (bCiteTng)
 226     {
 227         please_cite(fplog, "Lundborg2014");
 228     }
 229
 230     return of;
 231 }
 232
 233 ener_file_t mdoutf_get_fp_ene(gmx_mdoutf_t of)
 234 {
 235     return of->fp_ene;
 236 }
 237
 238 FILE *mdoutf_get_fp_dhdl(gmx_mdoutf_t of)
 239 {
 240     return of->fp_dhdl;
 241 }
 242
 243 gmx_wallcycle_t mdoutf_get_wcycle(gmx_mdoutf_t of)
 244 {
 245     return of->wcycle;
 246 }
 247
 248 void mdoutf_write_to_trajectory_files(FILE *fplog, const t_commrec *cr,
 249                                       gmx_mdoutf_t of,
 250                                       int mdof_flags,
 251                                       gmx_mtop_t *top_global,
 252                                       gmx_int64_t step, double t,
 253                                       t_state *state_local, t_state *state_global,
 254                                       ObservablesHistory *observablesHistory,
 255                                       gmx::ArrayRef<gmx::RVec> f_local)
 256 {
 257     rvec *f_global;
 258
 259     if (DOMAINDECOMP(cr))
 260     {
 261         if (mdof_flags & MDOF_CPT)
 262         {
 263             dd_collect_state(cr->dd, state_local, state_global);
 264         }
 265         else
 266         {
 267             if (mdof_flags & (MDOF_X | MDOF_X_COMPRESSED))
 268             {
 269                 gmx::ArrayRef<gmx::RVec> globalXRef = MASTER(cr) ? gmx::makeArrayRef(state_global->x) : gmx::EmptyArrayRef();
 270                 dd_collect_vec(cr->dd, state_local, state_local->x, globalXRef);
 271             }
 272             if (mdof_flags & MDOF_V)
 273             {
 274                 gmx::ArrayRef<gmx::RVec> globalVRef = MASTER(cr) ? gmx::makeArrayRef(state_global->v) : gmx::EmptyArrayRef();
 275                 dd_collect_vec(cr->dd, state_local, state_local->v, globalVRef);
 276             }
 277         }
 278         f_global = of->f_global;
 279         if (mdof_flags & MDOF_F)
 280         {
 281             dd_collect_vec(cr->dd, state_local, f_local, gmx::arrayRefFromArray(reinterpret_cast<gmx::RVec *>(f_global), f_local.size()));
 282         }
 283     }
 284     else
 285     {
 286         /* We have the whole state locally: copy the local state pointer */
 287         state_global = state_local;
 288
 289         f_global     = as_rvec_array(f_local.data());
 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             ivec one_ivec = { 1, 1, 1 };
 299             write_checkpoint(of->fn_cpt, of->bKeepAndNumCPT,
 300                              fplog, cr,
 301                              DOMAINDECOMP(cr) ? cr->dd->nc : one_ivec,
 302                              DOMAINDECOMP(cr) ? cr->dd->nnodes : cr->nnodes,
 303                              of->eIntegrator, of->simulation_part,
 304                              of->bExpanded, of->elamstats, step, t,
 305                              state_global, observablesHistory);
 306         }
 307
 308         if (mdof_flags & (MDOF_X | MDOF_V | MDOF_F))
 309         {
 310             const rvec *x = (mdof_flags & MDOF_X) ? as_rvec_array(state_global->x.data()) : nullptr;
 311             const rvec *v = (mdof_flags & MDOF_V) ? as_rvec_array(state_global->v.data()) : nullptr;
 312             const rvec *f = (mdof_flags & MDOF_F) ? f_global : nullptr;
 313
 314             if (of->fp_trn)
 315             {
 316                 gmx_trr_write_frame(of->fp_trn, step, t, state_local->lambda[efptFEP],
 317                                     state_local->box, top_global->natoms,
 318                                     x, v, f);
 319                 if (gmx_fio_flush(of->fp_trn) != 0)
 320                 {
 321                     gmx_file("Cannot write trajectory; maybe you are out of disk space?");
 322                 }
 323             }
 324
 325             /* If a TNG file is open for uncompressed coordinate output also write
 326                velocities and forces to it. */
 327             else if (of->tng)
 328             {
 329                 gmx_fwrite_tng(of->tng, FALSE, step, t, state_local->lambda[efptFEP],
 330                                state_local->box,
 331                                top_global->natoms,
 332                                x, v, f);
 333             }
 334             /* If only a TNG file is open for compressed coordinate output (no uncompressed
 335                coordinate output) also write forces and velocities to it. */
 336             else if (of->tng_low_prec)
 337             {
 338                 gmx_fwrite_tng(of->tng_low_prec, FALSE, step, t, state_local->lambda[efptFEP],
 339                                state_local->box,
 340                                top_global->natoms,
 341                                x, v, f);
 342             }
 343         }
 344         if (mdof_flags & MDOF_X_COMPRESSED)
 345         {
 346             rvec *xxtc = nullptr;
 347
 348             if (of->natoms_x_compressed == of->natoms_global)
 349             {
 350                 /* We are writing the positions of all of the atoms to
 351                    the compressed output */
 352                 xxtc = as_rvec_array(state_global->x.data());
 353             }
 354             else
 355             {
 356                 /* We are writing the positions of only a subset of
 357                    the atoms to the compressed output, so we have to
 358                    make a copy of the subset of coordinates. */
 359                 int i, j;
 360
 361                 snew(xxtc, of->natoms_x_compressed);
 362                 for (i = 0, j = 0; (i < of->natoms_global); i++)
 363                 {
 364                     if (ggrpnr(of->groups, egcCompressedX, i) == 0)
 365                     {
 366                         copy_rvec(state_global->x[i], xxtc[j++]);
 367                     }
 368                 }
 369             }
 370             if (write_xtc(of->fp_xtc, of->natoms_x_compressed, step, t,
 371                           state_local->box, xxtc, of->x_compression_precision) == 0)
 372             {
 373                 gmx_fatal(FARGS,
 374                           "XTC error. This indicates you are out of disk space, or a "
 375                           "simulation with major instabilities resulting in coordinates "
 376                           "that are NaN or too large to be represented in the XTC format.\n");
 377             }
 378             gmx_fwrite_tng(of->tng_low_prec,
 379                            TRUE,
 380                            step,
 381                            t,
 382                            state_local->lambda[efptFEP],
 383                            state_local->box,
 384                            of->natoms_x_compressed,
 385                            xxtc,
 386                            nullptr,
 387                            nullptr);
 388             if (of->natoms_x_compressed != of->natoms_global)
 389             {
 390                 sfree(xxtc);
 391             }
 392         }
 393         if (mdof_flags & (MDOF_BOX | MDOF_LAMBDA) && !(mdof_flags & (MDOF_X | MDOF_V | MDOF_F)) )
 394         {
 395             if (of->tng)
 396             {
 397                 real  lambda = -1;
 398                 rvec *box    = nullptr;
 399                 if (mdof_flags & MDOF_BOX)
 400                 {
 401                     box = state_local->box;
 402                 }
 403                 if (mdof_flags & MDOF_LAMBDA)
 404                 {
 405                     lambda = state_local->lambda[efptFEP];
 406                 }
 407                 gmx_fwrite_tng(of->tng, FALSE, step, t, lambda,
 408                                box, top_global->natoms,
 409                                nullptr, nullptr, nullptr);
 410             }
 411         }
 412         if (mdof_flags & (MDOF_BOX_COMPRESSED | MDOF_LAMBDA_COMPRESSED) && !(mdof_flags & (MDOF_X_COMPRESSED)) )
 413         {
 414             if (of->tng_low_prec)
 415             {
 416                 real  lambda = -1;
 417                 rvec *box    = nullptr;
 418                 if (mdof_flags & MDOF_BOX_COMPRESSED)
 419                 {
 420                     box = state_local->box;
 421                 }
 422                 if (mdof_flags & MDOF_LAMBDA_COMPRESSED)
 423                 {
 424                     lambda = state_local->lambda[efptFEP];
 425                 }
 426                 gmx_fwrite_tng(of->tng_low_prec, FALSE, step, t, lambda,
 427                                box, top_global->natoms,
 428                                nullptr, nullptr, nullptr);
 429             }
 430         }
 431     }
 432 }
 433
 434 void mdoutf_tng_close(gmx_mdoutf_t of)
 435 {
 436     if (of->tng || of->tng_low_prec)
 437     {
 438         wallcycle_start(of->wcycle, ewcTRAJ);
 439         gmx_tng_close(&of->tng);
 440         gmx_tng_close(&of->tng_low_prec);
 441         wallcycle_stop(of->wcycle, ewcTRAJ);
 442     }
 443 }
 444
 445 void done_mdoutf(gmx_mdoutf_t of)
 446 {
 447     if (of->fp_ene != nullptr)
 448     {
 449         close_enx(of->fp_ene);
 450     }
 451     if (of->fp_xtc)
 452     {
 453         close_xtc(of->fp_xtc);
 454     }
 455     if (of->fp_trn)
 456     {
 457         gmx_trr_close(of->fp_trn);
 458     }
 459     if (of->fp_dhdl != nullptr)
 460     {
 461         gmx_fio_fclose(of->fp_dhdl);
 462     }
 463     of->outputProvider->finishOutput();
 464     if (of->f_global != nullptr)
 465     {
 466         sfree(of->f_global);
 467     }
 468
 469     gmx_tng_close(&of->tng);
 470     gmx_tng_close(&of->tng_low_prec);
 471
 472     sfree(of);
 473 }
 474
 475 int mdoutf_get_tng_box_output_interval(gmx_mdoutf_t of)
 476 {
 477     if (of->tng)
 478     {
 479         return gmx_tng_get_box_output_interval(of->tng);
 480     }
 481     return 0;
 482 }
 483
 484 int mdoutf_get_tng_lambda_output_interval(gmx_mdoutf_t of)
 485 {
 486     if (of->tng)
 487     {
 488         return gmx_tng_get_lambda_output_interval(of->tng);
 489     }
 490     return 0;
 491 }
 492
 493 int mdoutf_get_tng_compressed_box_output_interval(gmx_mdoutf_t of)
 494 {
 495     if (of->tng_low_prec)
 496     {
 497         return gmx_tng_get_box_output_interval(of->tng_low_prec);
 498     }
 499     return 0;
 500 }
 501
 502 int mdoutf_get_tng_compressed_lambda_output_interval(gmx_mdoutf_t of)
 503 {
 504     if (of->tng_low_prec)
 505     {
 506         return gmx_tng_get_lambda_output_interval(of->tng_low_prec);
 507     }
 508     return 0;
 509 }