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