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 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.
   9  *
  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.
  14  *
  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.
  19  *
  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.
  24  *
  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.
  32  *
  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.
  35  */
  36 #include "gmxpre.h"
  37
  38 #include "mdoutf.h"
  39
  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"
  63
  64 struct gmx_mdoutf
  65 {
  66     t_fileio*                     fp_trn;
  67     t_fileio*                     fp_xtc;
  68     gmx_tng_trajectory_t          tng;
  69     gmx_tng_trajectory_t          tng_low_prec;
  70     int                           x_compression_precision; /* only used by XTC output */
  71     ener_file_t                   fp_ene;
  72     const char*                   fn_cpt;
  73     gmx_bool                      bKeepAndNumCPT;
  74     int                           eIntegrator;
  75     gmx_bool                      bExpanded;
  76     int                           elamstats;
  77     int                           simulation_part;
  78     FILE*                         fp_dhdl;
  79     int                           natoms_global;
  80     int                           natoms_x_compressed;
  81     SimulationGroups*             groups; /* for compressed position writing */
  82     gmx_wallcycle_t               wcycle;
  83     rvec*                         f_global;
  84     gmx::IMDOutputProvider*       outputProvider;
  85     const gmx::MdModulesNotifier* mdModulesNotifier;
  86 };
  87
  88
  89 gmx_mdoutf_t init_mdoutf(FILE*                         fplog,
  90                          int                           nfile,
  91                          const t_filenm                fnm[],
  92                          const gmx::MdrunOptions&      mdrunOptions,
  93                          const t_commrec*              cr,
  94                          gmx::IMDOutputProvider*       outputProvider,
  95                          const gmx::MdModulesNotifier& mdModulesNotifier,
  96                          const t_inputrec*             ir,
  97                          gmx_mtop_t*                   top_global,
  98                          const gmx_output_env_t*       oenv,
  99                          gmx_wallcycle_t               wcycle,
 100                          const gmx::StartingBehavior   startingBehavior)
 101 {
 102     gmx_mdoutf_t of;
 103     const char * appendMode = "a+", *writeMode = "w+", *filemode;
 104     gmx_bool     bCiteTng = FALSE;
 105     int          i;
 106     bool restartWithAppending = (startingBehavior == gmx::StartingBehavior::RestartWithAppending);
 107
 108     snew(of, 1);
 109
 110     of->fp_trn       = nullptr;
 111     of->fp_ene       = nullptr;
 112     of->fp_xtc       = nullptr;
 113     of->tng          = nullptr;
 114     of->tng_low_prec = nullptr;
 115     of->fp_dhdl      = nullptr;
 116
 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);
 122     of->wcycle                  = wcycle;
 123     of->f_global                = nullptr;
 124     of->outputProvider          = outputProvider;
 125
 126     if (MASTER(cr))
 127     {
 128         of->bKeepAndNumCPT = mdrunOptions.checkpointOptions.keepAndNumberCheckpointFiles;
 129
 130         filemode = restartWithAppending ? appendMode : writeMode;
 131
 132         if (EI_DYNAMICS(ir->eI) && ir->nstxout_compressed > 0)
 133         {
 134             const char* filename;
 135             filename = ftp2fn(efCOMPRESSED, nfile, fnm);
 136             switch (fn2ftp(filename))
 137             {
 138                 case efXTC: of->fp_xtc = open_xtc(filename, filemode); break;
 139                 case efTNG:
 140                     gmx_tng_open(filename, filemode[0], &of->tng_low_prec);
 141                     if (filemode[0] == 'w')
 142                     {
 143                         gmx_tng_prepare_low_prec_writing(of->tng_low_prec, top_global, ir);
 144                     }
 145                     bCiteTng = TRUE;
 146                     break;
 147                 default: gmx_incons("Invalid reduced precision file format");
 148             }
 149         }
 150         if ((EI_DYNAMICS(ir->eI) || EI_ENERGY_MINIMIZATION(ir->eI))
 151             && (!GMX_FAHCORE
 152                 && !(EI_DYNAMICS(ir->eI) && ir->nstxout == 0 && ir->nstvout == 0 && ir->nstfout == 0)))
 153         {
 154             const char* filename;
 155             filename = ftp2fn(efTRN, nfile, fnm);
 156             switch (fn2ftp(filename))
 157             {
 158                 case efTRR:
 159                 case efTRN:
 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)
 164                     {
 165                         of->fp_trn = gmx_trr_open(filename, filemode);
 166                     }
 167                     break;
 168                 case efTNG:
 169                     gmx_tng_open(filename, filemode[0], &of->tng);
 170                     if (filemode[0] == 'w')
 171                     {
 172                         gmx_tng_prepare_md_writing(of->tng, top_global, ir);
 173                     }
 174                     bCiteTng = TRUE;
 175                     break;
 176                 default: gmx_incons("Invalid full precision file format");
 177             }
 178         }
 179         if (EI_DYNAMICS(ir->eI) || EI_ENERGY_MINIMIZATION(ir->eI))
 180         {
 181             of->fp_ene = open_enx(ftp2fn(efEDR, nfile, fnm), filemode);
 182         }
 183         of->fn_cpt = opt2fn("-cpo", nfile, fnm);
 184
 185         if ((ir->efep != efepNO || ir->bSimTemp) && ir->fepvals->nstdhdl > 0
 186             && (ir->fepvals->separate_dhdl_file == esepdhdlfileYES) && EI_DYNAMICS(ir->eI))
 187         {
 188             if (restartWithAppending)
 189             {
 190                 of->fp_dhdl = gmx_fio_fopen(opt2fn("-dhdl", nfile, fnm), filemode);
 191             }
 192             else
 193             {
 194                 of->fp_dhdl = open_dhdl(opt2fn("-dhdl", nfile, fnm), ir, oenv);
 195             }
 196         }
 197
 198         outputProvider->initOutput(fplog, nfile, fnm, restartWithAppending, oenv);
 199         of->mdModulesNotifier = &mdModulesNotifier;
 200
 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++)
 209         {
 210             if (getGroupType(*of->groups, SimulationAtomGroupType::CompressedPositionOutput, i) == 0)
 211             {
 212                 of->natoms_x_compressed++;
 213             }
 214         }
 215
 216         if (ir->nstfout && DOMAINDECOMP(cr))
 217         {
 218             snew(of->f_global, top_global->natoms);
 219         }
 220     }
 221
 222     if (bCiteTng)
 223     {
 224         please_cite(fplog, "Lundborg2014");
 225     }
 226
 227     return of;
 228 }
 229
 230 ener_file_t mdoutf_get_fp_ene(gmx_mdoutf_t of)
 231 {
 232     return of->fp_ene;
 233 }
 234
 235 FILE* mdoutf_get_fp_dhdl(gmx_mdoutf_t of)
 236 {
 237     return of->fp_dhdl;
 238 }
 239
 240 gmx_wallcycle_t mdoutf_get_wcycle(gmx_mdoutf_t of)
 241 {
 242     return of->wcycle;
 243 }
 244
 245 void mdoutf_write_to_trajectory_files(FILE*                    fplog,
 246                                       const t_commrec*         cr,
 247                                       gmx_mdoutf_t             of,
 248                                       int                      mdof_flags,
 249                                       int                      natoms,
 250                                       int64_t                  step,
 251                                       double                   t,
 252                                       t_state*                 state_local,
 253                                       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                 auto globalXRef = MASTER(cr) ? state_global->x : gmx::ArrayRef<gmx::RVec>();
 270                 dd_collect_vec(cr->dd, state_local, state_local->x, globalXRef);
 271             }
 272             if (mdof_flags & MDOF_V)
 273             {
 274                 auto globalVRef = MASTER(cr) ? state_global->v : gmx::ArrayRef<gmx::RVec>();
 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,
 282                            gmx::arrayRefFromArray(reinterpret_cast<gmx::RVec*>(f_global), f_local.size()));
 283         }
 284     }
 285     else
 286     {
 287         /* We have the whole state locally: copy the local state pointer */
 288         state_global = state_local;
 289
 290         f_global = as_rvec_array(f_local.data());
 291     }
 292
 293     if (MASTER(cr))
 294     {
 295         if (mdof_flags & MDOF_CPT)
 296         {
 297             fflush_tng(of->tng);
 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));
 305         }
 306
 307         if (mdof_flags & (MDOF_X | MDOF_V | MDOF_F))
 308         {
 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;
 312
 313             if (of->fp_trn)
 314             {
 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)
 318                 {
 319                     gmx_file("Cannot write trajectory; maybe you are out of disk space?");
 320                 }
 321             }
 322
 323             /* If a TNG file is open for uncompressed coordinate output also write
 324                velocities and forces to it. */
 325             else if (of->tng)
 326             {
 327                 gmx_fwrite_tng(of->tng, FALSE, step, t, state_local->lambda[efptFEP],
 328                                state_local->box, natoms, x, v, f);
 329             }
 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)
 333             {
 334                 gmx_fwrite_tng(of->tng_low_prec, FALSE, step, t, state_local->lambda[efptFEP],
 335                                state_local->box, natoms, x, v, f);
 336             }
 337         }
 338         if (mdof_flags & MDOF_X_COMPRESSED)
 339         {
 340             rvec* xxtc = nullptr;
 341
 342             if (of->natoms_x_compressed == of->natoms_global)
 343             {
 344                 /* We are writing the positions of all of the atoms to
 345                    the compressed output */
 346                 xxtc = state_global->x.rvec_array();
 347             }
 348             else
 349             {
 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. */
 353                 int i, j;
 354
 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++)
 358                 {
 359                     if (getGroupType(*of->groups, SimulationAtomGroupType::CompressedPositionOutput, i) == 0)
 360                     {
 361                         copy_rvec(x[i], xxtc[j++]);
 362                     }
 363                 }
 364             }
 365             if (write_xtc(of->fp_xtc, of->natoms_x_compressed, step, t, state_local->box, xxtc,
 366                           of->x_compression_precision)
 367                 == 0)
 368             {
 369                 gmx_fatal(FARGS,
 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");
 373             }
 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)
 377             {
 378                 sfree(xxtc);
 379             }
 380         }
 381         if (mdof_flags & (MDOF_BOX | MDOF_LAMBDA) && !(mdof_flags & (MDOF_X | MDOF_V | MDOF_F)))
 382         {
 383             if (of->tng)
 384             {
 385                 real  lambda = -1;
 386                 rvec* box    = nullptr;
 387                 if (mdof_flags & MDOF_BOX)
 388                 {
 389                     box = state_local->box;
 390                 }
 391                 if (mdof_flags & MDOF_LAMBDA)
 392                 {
 393                     lambda = state_local->lambda[efptFEP];
 394                 }
 395                 gmx_fwrite_tng(of->tng, FALSE, step, t, lambda, box, natoms, nullptr, nullptr, nullptr);
 396             }
 397         }
 398         if (mdof_flags & (MDOF_BOX_COMPRESSED | MDOF_LAMBDA_COMPRESSED)
 399             && !(mdof_flags & (MDOF_X_COMPRESSED)))
 400         {
 401             if (of->tng_low_prec)
 402             {
 403                 real  lambda = -1;
 404                 rvec* box    = nullptr;
 405                 if (mdof_flags & MDOF_BOX_COMPRESSED)
 406                 {
 407                     box = state_local->box;
 408                 }
 409                 if (mdof_flags & MDOF_LAMBDA_COMPRESSED)
 410                 {
 411                     lambda = state_local->lambda[efptFEP];
 412                 }
 413                 gmx_fwrite_tng(of->tng_low_prec, FALSE, step, t, lambda, box, natoms, nullptr,
 414                                nullptr, nullptr);
 415             }
 416         }
 417     }
 418 }
 419
 420 void mdoutf_tng_close(gmx_mdoutf_t of)
 421 {
 422     if (of->tng || of->tng_low_prec)
 423     {
 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);
 428     }
 429 }
 430
 431 void done_mdoutf(gmx_mdoutf_t of)
 432 {
 433     if (of->fp_ene != nullptr)
 434     {
 435         done_ener_file(of->fp_ene);
 436     }
 437     if (of->fp_xtc)
 438     {
 439         close_xtc(of->fp_xtc);
 440     }
 441     if (of->fp_trn)
 442     {
 443         gmx_trr_close(of->fp_trn);
 444     }
 445     if (of->fp_dhdl != nullptr)
 446     {
 447         gmx_fio_fclose(of->fp_dhdl);
 448     }
 449     of->outputProvider->finishOutput();
 450     if (of->f_global != nullptr)
 451     {
 452         sfree(of->f_global);
 453     }
 454
 455     gmx_tng_close(&of->tng);
 456     gmx_tng_close(&of->tng_low_prec);
 457
 458     sfree(of);
 459 }
 460
 461 int mdoutf_get_tng_box_output_interval(gmx_mdoutf_t of)
 462 {
 463     if (of->tng)
 464     {
 465         return gmx_tng_get_box_output_interval(of->tng);
 466     }
 467     return 0;
 468 }
 469
 470 int mdoutf_get_tng_lambda_output_interval(gmx_mdoutf_t of)
 471 {
 472     if (of->tng)
 473     {
 474         return gmx_tng_get_lambda_output_interval(of->tng);
 475     }
 476     return 0;
 477 }
 478
 479 int mdoutf_get_tng_compressed_box_output_interval(gmx_mdoutf_t of)
 480 {
 481     if (of->tng_low_prec)
 482     {
 483         return gmx_tng_get_box_output_interval(of->tng_low_prec);
 484     }
 485     return 0;
 486 }
 487
 488 int mdoutf_get_tng_compressed_lambda_output_interval(gmx_mdoutf_t of)
 489 {
 490     if (of->tng_low_prec)
 491     {
 492         return gmx_tng_get_lambda_output_interval(of->tng_low_prec);
 493     }
 494     return 0;
 495 }