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,2021, 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 "config.h"
  41
  42 #include "gromacs/commandline/filenm.h"
  43 #include "gromacs/domdec/collect.h"
  44 #include "gromacs/domdec/domdec_struct.h"
  45 #include "gromacs/fileio/checkpoint.h"
  46 #include "gromacs/fileio/gmxfio.h"
  47 #include "gromacs/fileio/tngio.h"
  48 #include "gromacs/fileio/trrio.h"
  49 #include "gromacs/fileio/xtcio.h"
  50 #include "gromacs/fileio/xvgr.h"
  51 #include "gromacs/math/vec.h"
  52 #include "gromacs/mdlib/trajectory_writing.h"
  53 #include "gromacs/mdrunutility/handlerestart.h"
  54 #include "gromacs/mdrunutility/multisim.h"
  55 #include "gromacs/mdtypes/awh_history.h"
  56 #include "gromacs/mdtypes/commrec.h"
  57 #include "gromacs/mdtypes/df_history.h"
  58 #include "gromacs/mdtypes/edsamhistory.h"
  59 #include "gromacs/mdtypes/energyhistory.h"
  60 #include "gromacs/mdtypes/imdoutputprovider.h"
  61 #include "gromacs/mdtypes/inputrec.h"
  62 #include "gromacs/mdtypes/md_enums.h"
  63 #include "gromacs/mdtypes/mdrunoptions.h"
  64 #include "gromacs/mdtypes/observableshistory.h"
  65 #include "gromacs/mdtypes/state.h"
  66 #include "gromacs/mdtypes/swaphistory.h"
  67 #include "gromacs/timing/wallcycle.h"
  68 #include "gromacs/topology/topology.h"
  69 #include "gromacs/utility/baseversion.h"
  70 #include "gromacs/utility/fatalerror.h"
  71 #include "gromacs/utility/pleasecite.h"
  72 #include "gromacs/utility/programcontext.h"
  73 #include "gromacs/utility/smalloc.h"
  74 #include "gromacs/utility/sysinfo.h"
  75
  76 struct gmx_mdoutf
  77 {
  78     t_fileio*                     fp_trn;
  79     t_fileio*                     fp_xtc;
  80     gmx_tng_trajectory_t          tng;
  81     gmx_tng_trajectory_t          tng_low_prec;
  82     int                           x_compression_precision; /* only used by XTC output */
  83     ener_file_t                   fp_ene;
  84     const char*                   fn_cpt;
  85     gmx_bool                      bKeepAndNumCPT;
  86     IntegrationAlgorithm          eIntegrator;
  87     gmx_bool                      bExpanded;
  88     LambdaWeightCalculation       elamstats;
  89     int                           simulation_part;
  90     FILE*                         fp_dhdl;
  91     int                           natoms_global;
  92     int                           natoms_x_compressed;
  93     const SimulationGroups*       groups; /* for compressed position writing */
  94     gmx_wallcycle_t               wcycle;
  95     rvec*                         f_global;
  96     gmx::IMDOutputProvider*       outputProvider;
  97     const gmx::MdModulesNotifier* mdModulesNotifier;
  98     bool                          simulationsShareState;
  99     MPI_Comm                      mastersComm;
 100 };
 101
 102
 103 gmx_mdoutf_t init_mdoutf(FILE*                         fplog,
 104                          int                           nfile,
 105                          const t_filenm                fnm[],
 106                          const gmx::MdrunOptions&      mdrunOptions,
 107                          const t_commrec*              cr,
 108                          gmx::IMDOutputProvider*       outputProvider,
 109                          const gmx::MdModulesNotifier& mdModulesNotifier,
 110                          const t_inputrec*             ir,
 111                          const gmx_mtop_t*             top_global,
 112                          const gmx_output_env_t*       oenv,
 113                          gmx_wallcycle_t               wcycle,
 114                          const gmx::StartingBehavior   startingBehavior,
 115                          bool                          simulationsShareState,
 116                          const gmx_multisim_t*         ms)
 117 {
 118     gmx_mdoutf_t of;
 119     const char * appendMode = "a+", *writeMode = "w+", *filemode;
 120     gmx_bool     bCiteTng = FALSE;
 121     int          i;
 122     bool restartWithAppending = (startingBehavior == gmx::StartingBehavior::RestartWithAppending);
 123
 124     snew(of, 1);
 125
 126     of->fp_trn       = nullptr;
 127     of->fp_ene       = nullptr;
 128     of->fp_xtc       = nullptr;
 129     of->tng          = nullptr;
 130     of->tng_low_prec = nullptr;
 131     of->fp_dhdl      = nullptr;
 132
 133     of->eIntegrator             = ir->eI;
 134     of->bExpanded               = ir->bExpanded;
 135     of->elamstats               = ir->expandedvals->elamstats;
 136     of->simulation_part         = ir->simulation_part;
 137     of->x_compression_precision = static_cast<int>(ir->x_compression_precision);
 138     of->wcycle                  = wcycle;
 139     of->f_global                = nullptr;
 140     of->outputProvider          = outputProvider;
 141
 142     GMX_RELEASE_ASSERT(!simulationsShareState || ms != nullptr,
 143                        "Need valid multisim object when simulations share state");
 144     of->simulationsShareState = simulationsShareState;
 145     if (of->simulationsShareState)
 146     {
 147         of->mastersComm = ms->mastersComm_;
 148     }
 149
 150     if (MASTER(cr))
 151     {
 152         of->bKeepAndNumCPT = mdrunOptions.checkpointOptions.keepAndNumberCheckpointFiles;
 153
 154         filemode = restartWithAppending ? appendMode : writeMode;
 155
 156         if (EI_DYNAMICS(ir->eI) && ir->nstxout_compressed > 0)
 157         {
 158             const char* filename;
 159             filename = ftp2fn(efCOMPRESSED, nfile, fnm);
 160             switch (fn2ftp(filename))
 161             {
 162                 case efXTC: of->fp_xtc = open_xtc(filename, filemode); break;
 163                 case efTNG:
 164                     gmx_tng_open(filename, filemode[0], &of->tng_low_prec);
 165                     if (filemode[0] == 'w')
 166                     {
 167                         gmx_tng_prepare_low_prec_writing(of->tng_low_prec, top_global, ir);
 168                     }
 169                     bCiteTng = TRUE;
 170                     break;
 171                 default: gmx_incons("Invalid reduced precision file format");
 172             }
 173         }
 174         if ((EI_DYNAMICS(ir->eI) || EI_ENERGY_MINIMIZATION(ir->eI))
 175             && (!GMX_FAHCORE
 176                 && !(EI_DYNAMICS(ir->eI) && ir->nstxout == 0 && ir->nstvout == 0 && ir->nstfout == 0)))
 177         {
 178             const char* filename;
 179             filename = ftp2fn(efTRN, nfile, fnm);
 180             switch (fn2ftp(filename))
 181             {
 182                 case efTRR:
 183                 case efTRN:
 184                     /* If there is no uncompressed coordinate output and
 185                        there is compressed TNG output write forces
 186                        and/or velocities to the TNG file instead. */
 187                     if (ir->nstxout != 0 || ir->nstxout_compressed == 0 || !of->tng_low_prec)
 188                     {
 189                         of->fp_trn = gmx_trr_open(filename, filemode);
 190                     }
 191                     break;
 192                 case efTNG:
 193                     gmx_tng_open(filename, filemode[0], &of->tng);
 194                     if (filemode[0] == 'w')
 195                     {
 196                         gmx_tng_prepare_md_writing(of->tng, top_global, ir);
 197                     }
 198                     bCiteTng = TRUE;
 199                     break;
 200                 default: gmx_incons("Invalid full precision file format");
 201             }
 202         }
 203         if (EI_DYNAMICS(ir->eI) || EI_ENERGY_MINIMIZATION(ir->eI))
 204         {
 205             of->fp_ene = open_enx(ftp2fn(efEDR, nfile, fnm), filemode);
 206         }
 207         of->fn_cpt = opt2fn("-cpo", nfile, fnm);
 208
 209         if ((ir->efep != FreeEnergyPerturbationType::No || ir->bSimTemp) && ir->fepvals->nstdhdl > 0
 210             && (ir->fepvals->separate_dhdl_file == SeparateDhdlFile::Yes) && EI_DYNAMICS(ir->eI))
 211         {
 212             if (restartWithAppending)
 213             {
 214                 of->fp_dhdl = gmx_fio_fopen(opt2fn("-dhdl", nfile, fnm), filemode);
 215             }
 216             else
 217             {
 218                 of->fp_dhdl = open_dhdl(opt2fn("-dhdl", nfile, fnm), ir, oenv);
 219             }
 220         }
 221
 222         outputProvider->initOutput(fplog, nfile, fnm, restartWithAppending, oenv);
 223         of->mdModulesNotifier = &mdModulesNotifier;
 224
 225         /* Set up atom counts so they can be passed to actual
 226            trajectory-writing routines later. Also, XTC writing needs
 227            to know what (and how many) atoms might be in the XTC
 228            groups, and how to look up later which ones they are. */
 229         of->natoms_global       = top_global->natoms;
 230         of->groups              = &top_global->groups;
 231         of->natoms_x_compressed = 0;
 232         for (i = 0; (i < top_global->natoms); i++)
 233         {
 234             if (getGroupType(*of->groups, SimulationAtomGroupType::CompressedPositionOutput, i) == 0)
 235             {
 236                 of->natoms_x_compressed++;
 237             }
 238         }
 239
 240         if (ir->nstfout && DOMAINDECOMP(cr))
 241         {
 242             snew(of->f_global, top_global->natoms);
 243         }
 244     }
 245
 246     if (bCiteTng)
 247     {
 248         please_cite(fplog, "Lundborg2014");
 249     }
 250
 251     return of;
 252 }
 253
 254 ener_file_t mdoutf_get_fp_ene(gmx_mdoutf_t of)
 255 {
 256     return of->fp_ene;
 257 }
 258
 259 FILE* mdoutf_get_fp_dhdl(gmx_mdoutf_t of)
 260 {
 261     return of->fp_dhdl;
 262 }
 263
 264 gmx_wallcycle_t mdoutf_get_wcycle(gmx_mdoutf_t of)
 265 {
 266     return of->wcycle;
 267 }
 268
 269 static void mpiBarrierBeforeRename(const bool applyMpiBarrierBeforeRename, MPI_Comm mpiBarrierCommunicator)
 270 {
 271     if (applyMpiBarrierBeforeRename)
 272     {
 273 #if GMX_MPI
 274         MPI_Barrier(mpiBarrierCommunicator);
 275 #else
 276         GMX_RELEASE_ASSERT(false, "Should not request a barrier without MPI");
 277         GMX_UNUSED_VALUE(mpiBarrierCommunicator);
 278 #endif
 279     }
 280 }
 281 /*! \brief Write a checkpoint to the filename
 282  *
 283  * Appends the _step<step>.cpt with bNumberAndKeep, otherwise moves
 284  * the previous checkpoint filename with suffix _prev.cpt.
 285  */
 286 static void write_checkpoint(const char*                     fn,
 287                              gmx_bool                        bNumberAndKeep,
 288                              FILE*                           fplog,
 289                              const t_commrec*                cr,
 290                              ivec                            domdecCells,
 291                              int                             nppnodes,
 292                              IntegrationAlgorithm            eIntegrator,
 293                              int                             simulation_part,
 294                              gmx_bool                        bExpanded,
 295                              LambdaWeightCalculation         elamstats,
 296                              int64_t                         step,
 297                              double                          t,
 298                              t_state*                        state,
 299                              ObservablesHistory*             observablesHistory,
 300                              const gmx::MdModulesNotifier&   mdModulesNotifier,
 301                              gmx::WriteCheckpointDataHolder* modularSimulatorCheckpointData,
 302                              bool                            applyMpiBarrierBeforeRename,
 303                              MPI_Comm                        mpiBarrierCommunicator)
 304 {
 305     t_fileio* fp;
 306     char*     fntemp; /* the temporary checkpoint file name */
 307     int       npmenodes;
 308     char      buf[1024], suffix[5 + STEPSTRSIZE], sbuf[STEPSTRSIZE];
 309     t_fileio* ret;
 310
 311     if (DOMAINDECOMP(cr))
 312     {
 313         npmenodes = cr->npmenodes;
 314     }
 315     else
 316     {
 317         npmenodes = 0;
 318     }
 319
 320 #if !GMX_NO_RENAME
 321     /* make the new temporary filename */
 322     snew(fntemp, std::strlen(fn) + 5 + STEPSTRSIZE);
 323     std::strcpy(fntemp, fn);
 324     fntemp[std::strlen(fn) - std::strlen(ftp2ext(fn2ftp(fn))) - 1] = '\0';
 325     sprintf(suffix, "_%s%s", "step", gmx_step_str(step, sbuf));
 326     std::strcat(fntemp, suffix);
 327     std::strcat(fntemp, fn + std::strlen(fn) - std::strlen(ftp2ext(fn2ftp(fn))) - 1);
 328 #else
 329     /* if we can't rename, we just overwrite the cpt file.
 330      * dangerous if interrupted.
 331      */
 332     snew(fntemp, std::strlen(fn));
 333     std::strcpy(fntemp, fn);
 334 #endif
 335     std::string timebuf = gmx_format_current_time();
 336
 337     if (fplog)
 338     {
 339         fprintf(fplog, "Writing checkpoint, step %s at %s\n\n", gmx_step_str(step, buf), timebuf.c_str());
 340     }
 341
 342     /* Get offsets for open files */
 343     auto outputfiles = gmx_fio_get_output_file_positions();
 344
 345     fp = gmx_fio_open(fntemp, "w");
 346
 347     /* We can check many more things now (CPU, acceleration, etc), but
 348      * it is highly unlikely to have two separate builds with exactly
 349      * the same version, user, time, and build host!
 350      */
 351
 352     int nlambda = (state->dfhist ? state->dfhist->nlambda : 0);
 353
 354     edsamhistory_t* edsamhist = observablesHistory->edsamHistory.get();
 355     int             nED       = (edsamhist ? edsamhist->nED : 0);
 356
 357     swaphistory_t* swaphist    = observablesHistory->swapHistory.get();
 358     SwapType       eSwapCoords = (swaphist ? swaphist->eSwapCoords : SwapType::No);
 359
 360     CheckpointHeaderContents headerContents = { 0,
 361                                                 { 0 },
 362                                                 { 0 },
 363                                                 { 0 },
 364                                                 { 0 },
 365                                                 GMX_DOUBLE,
 366                                                 { 0 },
 367                                                 { 0 },
 368                                                 eIntegrator,
 369                                                 simulation_part,
 370                                                 step,
 371                                                 t,
 372                                                 nppnodes,
 373                                                 { 0 },
 374                                                 npmenodes,
 375                                                 state->natoms,
 376                                                 state->ngtc,
 377                                                 state->nnhpres,
 378                                                 state->nhchainlength,
 379                                                 nlambda,
 380                                                 state->flags,
 381                                                 0,
 382                                                 0,
 383                                                 0,
 384                                                 0,
 385                                                 0,
 386                                                 nED,
 387                                                 eSwapCoords,
 388                                                 false };
 389     std::strcpy(headerContents.version, gmx_version());
 390     std::strcpy(headerContents.fprog, gmx::getProgramContext().fullBinaryPath());
 391     std::strcpy(headerContents.ftime, timebuf.c_str());
 392     if (DOMAINDECOMP(cr))
 393     {
 394         copy_ivec(domdecCells, headerContents.dd_nc);
 395     }
 396
 397     write_checkpoint_data(fp,
 398                           headerContents,
 399                           bExpanded,
 400                           elamstats,
 401                           state,
 402                           observablesHistory,
 403                           mdModulesNotifier,
 404                           &outputfiles,
 405                           modularSimulatorCheckpointData);
 406
 407     /* we really, REALLY, want to make sure to physically write the checkpoint,
 408        and all the files it depends on, out to disk. Because we've
 409        opened the checkpoint with gmx_fio_open(), it's in our list
 410        of open files.  */
 411     ret = gmx_fio_all_output_fsync();
 412
 413     if (ret)
 414     {
 415         char buf[STRLEN];
 416         sprintf(buf, "Cannot fsync '%s'; maybe you are out of disk space?", gmx_fio_getname(ret));
 417
 418         if (getenv(GMX_IGNORE_FSYNC_FAILURE_ENV) == nullptr)
 419         {
 420             gmx_file(buf);
 421         }
 422         else
 423         {
 424             gmx_warning("%s", buf);
 425         }
 426     }
 427
 428     if (gmx_fio_close(fp) != 0)
 429     {
 430         gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");
 431     }
 432
 433     /* we don't move the checkpoint if the user specified they didn't want it,
 434        or if the fsyncs failed */
 435 #if !GMX_NO_RENAME
 436     if (!bNumberAndKeep && !ret)
 437     {
 438         if (gmx_fexist(fn))
 439         {
 440             /* Rename the previous checkpoint file */
 441             mpiBarrierBeforeRename(applyMpiBarrierBeforeRename, mpiBarrierCommunicator);
 442
 443             std::strcpy(buf, fn);
 444             buf[std::strlen(fn) - std::strlen(ftp2ext(fn2ftp(fn))) - 1] = '\0';
 445             std::strcat(buf, "_prev");
 446             std::strcat(buf, fn + std::strlen(fn) - std::strlen(ftp2ext(fn2ftp(fn))) - 1);
 447             if (!GMX_FAHCORE)
 448             {
 449                 /* we copy here so that if something goes wrong between now and
 450                  * the rename below, there's always a state.cpt.
 451                  * If renames are atomic (such as in POSIX systems),
 452                  * this copying should be unneccesary.
 453                  */
 454                 gmx_file_copy(fn, buf, FALSE);
 455                 /* We don't really care if this fails:
 456                  * there's already a new checkpoint.
 457                  */
 458             }
 459             else
 460             {
 461                 gmx_file_rename(fn, buf);
 462             }
 463         }
 464
 465         /* Rename the checkpoint file from the temporary to the final name */
 466         mpiBarrierBeforeRename(applyMpiBarrierBeforeRename, mpiBarrierCommunicator);
 467
 468         if (gmx_file_rename(fntemp, fn) != 0)
 469         {
 470             gmx_file("Cannot rename checkpoint file; maybe you are out of disk space?");
 471         }
 472     }
 473 #endif /* GMX_NO_RENAME */
 474
 475     sfree(fntemp);
 476
 477 #if GMX_FAHCORE
 478     /*code for alternate checkpointing scheme.  moved from top of loop over
 479        steps */
 480     fcRequestCheckPoint();
 481     if (fcCheckPointParallel(cr->nodeid, NULL, 0) == 0)
 482     {
 483         gmx_fatal(3, __FILE__, __LINE__, "Checkpoint error on step %d\n", step);
 484     }
 485 #endif /* end GMX_FAHCORE block */
 486 }
 487
 488 void mdoutf_write_checkpoint(gmx_mdoutf_t                    of,
 489                              FILE*                           fplog,
 490                              const t_commrec*                cr,
 491                              int64_t                         step,
 492                              double                          t,
 493                              t_state*                        state_global,
 494                              ObservablesHistory*             observablesHistory,
 495                              gmx::WriteCheckpointDataHolder* modularSimulatorCheckpointData)
 496 {
 497     fflush_tng(of->tng);
 498     fflush_tng(of->tng_low_prec);
 499     /* Write the checkpoint file.
 500      * When simulations share the state, an MPI barrier is applied before
 501      * renaming old and new checkpoint files to minimize the risk of
 502      * checkpoint files getting out of sync.
 503      */
 504     ivec one_ivec = { 1, 1, 1 };
 505     write_checkpoint(of->fn_cpt,
 506                      of->bKeepAndNumCPT,
 507                      fplog,
 508                      cr,
 509                      DOMAINDECOMP(cr) ? cr->dd->numCells : one_ivec,
 510                      DOMAINDECOMP(cr) ? cr->dd->nnodes : cr->nnodes,
 511                      of->eIntegrator,
 512                      of->simulation_part,
 513                      of->bExpanded,
 514                      of->elamstats,
 515                      step,
 516                      t,
 517                      state_global,
 518                      observablesHistory,
 519                      *(of->mdModulesNotifier),
 520                      modularSimulatorCheckpointData,
 521                      of->simulationsShareState,
 522                      of->mastersComm);
 523 }
 524
 525 void mdoutf_write_to_trajectory_files(FILE*                           fplog,
 526                                       const t_commrec*                cr,
 527                                       gmx_mdoutf_t                    of,
 528                                       int                             mdof_flags,
 529                                       int                             natoms,
 530                                       int64_t                         step,
 531                                       double                          t,
 532                                       t_state*                        state_local,
 533                                       t_state*                        state_global,
 534                                       ObservablesHistory*             observablesHistory,
 535                                       gmx::ArrayRef<const gmx::RVec>  f_local,
 536                                       gmx::WriteCheckpointDataHolder* modularSimulatorCheckpointData)
 537 {
 538     const rvec* f_global;
 539
 540     if (DOMAINDECOMP(cr))
 541     {
 542         if (mdof_flags & MDOF_CPT)
 543         {
 544             dd_collect_state(cr->dd, state_local, state_global);
 545         }
 546         else
 547         {
 548             if (mdof_flags & (MDOF_X | MDOF_X_COMPRESSED))
 549             {
 550                 auto globalXRef = MASTER(cr) ? state_global->x : gmx::ArrayRef<gmx::RVec>();
 551                 dd_collect_vec(cr->dd,
 552                                state_local->ddp_count,
 553                                state_local->ddp_count_cg_gl,
 554                                state_local->cg_gl,
 555                                state_local->x,
 556                                globalXRef);
 557             }
 558             if (mdof_flags & MDOF_V)
 559             {
 560                 auto globalVRef = MASTER(cr) ? state_global->v : gmx::ArrayRef<gmx::RVec>();
 561                 dd_collect_vec(cr->dd,
 562                                state_local->ddp_count,
 563                                state_local->ddp_count_cg_gl,
 564                                state_local->cg_gl,
 565                                state_local->v,
 566                                globalVRef);
 567             }
 568         }
 569         f_global = of->f_global;
 570         if (mdof_flags & MDOF_F)
 571         {
 572             auto globalFRef =
 573                     MASTER(cr) ? gmx::arrayRefFromArray(reinterpret_cast<gmx::RVec*>(of->f_global),
 574                                                         f_local.size())
 575                                : gmx::ArrayRef<gmx::RVec>();
 576             dd_collect_vec(cr->dd,
 577                            state_local->ddp_count,
 578                            state_local->ddp_count_cg_gl,
 579                            state_local->cg_gl,
 580                            f_local,
 581                            globalFRef);
 582         }
 583     }
 584     else
 585     {
 586         /* We have the whole state locally: copy the local state pointer */
 587         state_global = state_local;
 588
 589         f_global = as_rvec_array(f_local.data());
 590     }
 591
 592     if (MASTER(cr))
 593     {
 594         if (mdof_flags & MDOF_CPT)
 595         {
 596             mdoutf_write_checkpoint(
 597                     of, fplog, cr, step, t, state_global, observablesHistory, modularSimulatorCheckpointData);
 598         }
 599
 600         if (mdof_flags & (MDOF_X | MDOF_V | MDOF_F))
 601         {
 602             const rvec* x = (mdof_flags & MDOF_X) ? state_global->x.rvec_array() : nullptr;
 603             const rvec* v = (mdof_flags & MDOF_V) ? state_global->v.rvec_array() : nullptr;
 604             const rvec* f = (mdof_flags & MDOF_F) ? f_global : nullptr;
 605
 606             if (of->fp_trn)
 607             {
 608                 gmx_trr_write_frame(of->fp_trn,
 609                                     step,
 610                                     t,
 611                                     state_local->lambda[FreeEnergyPerturbationCouplingType::Fep],
 612                                     state_local->box,
 613                                     natoms,
 614                                     x,
 615                                     v,
 616                                     f);
 617                 if (gmx_fio_flush(of->fp_trn) != 0)
 618                 {
 619                     gmx_file("Cannot write trajectory; maybe you are out of disk space?");
 620                 }
 621             }
 622
 623             /* If a TNG file is open for uncompressed coordinate output also write
 624                velocities and forces to it. */
 625             else if (of->tng)
 626             {
 627                 gmx_fwrite_tng(of->tng,
 628                                FALSE,
 629                                step,
 630                                t,
 631                                state_local->lambda[FreeEnergyPerturbationCouplingType::Fep],
 632                                state_local->box,
 633                                natoms,
 634                                x,
 635                                v,
 636                                f);
 637             }
 638             /* If only a TNG file is open for compressed coordinate output (no uncompressed
 639                coordinate output) also write forces and velocities to it. */
 640             else if (of->tng_low_prec)
 641             {
 642                 gmx_fwrite_tng(of->tng_low_prec,
 643                                FALSE,
 644                                step,
 645                                t,
 646                                state_local->lambda[FreeEnergyPerturbationCouplingType::Fep],
 647                                state_local->box,
 648                                natoms,
 649                                x,
 650                                v,
 651                                f);
 652             }
 653         }
 654         if (mdof_flags & MDOF_X_COMPRESSED)
 655         {
 656             rvec* xxtc = nullptr;
 657
 658             if (of->natoms_x_compressed == of->natoms_global)
 659             {
 660                 /* We are writing the positions of all of the atoms to
 661                    the compressed output */
 662                 xxtc = state_global->x.rvec_array();
 663             }
 664             else
 665             {
 666                 /* We are writing the positions of only a subset of
 667                    the atoms to the compressed output, so we have to
 668                    make a copy of the subset of coordinates. */
 669                 int i, j;
 670
 671                 snew(xxtc, of->natoms_x_compressed);
 672                 auto x = makeArrayRef(state_global->x);
 673                 for (i = 0, j = 0; (i < of->natoms_global); i++)
 674                 {
 675                     if (getGroupType(*of->groups, SimulationAtomGroupType::CompressedPositionOutput, i) == 0)
 676                     {
 677                         copy_rvec(x[i], xxtc[j++]);
 678                     }
 679                 }
 680             }
 681             if (write_xtc(of->fp_xtc, of->natoms_x_compressed, step, t, state_local->box, xxtc, of->x_compression_precision)
 682                 == 0)
 683             {
 684                 gmx_fatal(FARGS,
 685                           "XTC error. This indicates you are out of disk space, or a "
 686                           "simulation with major instabilities resulting in coordinates "
 687                           "that are NaN or too large to be represented in the XTC format.\n");
 688             }
 689             gmx_fwrite_tng(of->tng_low_prec,
 690                            TRUE,
 691                            step,
 692                            t,
 693                            state_local->lambda[FreeEnergyPerturbationCouplingType::Fep],
 694                            state_local->box,
 695                            of->natoms_x_compressed,
 696                            xxtc,
 697                            nullptr,
 698                            nullptr);
 699             if (of->natoms_x_compressed != of->natoms_global)
 700             {
 701                 sfree(xxtc);
 702             }
 703         }
 704         if (mdof_flags & (MDOF_BOX | MDOF_LAMBDA) && !(mdof_flags & (MDOF_X | MDOF_V | MDOF_F)))
 705         {
 706             if (of->tng)
 707             {
 708                 real  lambda = -1;
 709                 rvec* box    = nullptr;
 710                 if (mdof_flags & MDOF_BOX)
 711                 {
 712                     box = state_local->box;
 713                 }
 714                 if (mdof_flags & MDOF_LAMBDA)
 715                 {
 716                     lambda = state_local->lambda[FreeEnergyPerturbationCouplingType::Fep];
 717                 }
 718                 gmx_fwrite_tng(of->tng, FALSE, step, t, lambda, box, natoms, nullptr, nullptr, nullptr);
 719             }
 720         }
 721         if (mdof_flags & (MDOF_BOX_COMPRESSED | MDOF_LAMBDA_COMPRESSED)
 722             && !(mdof_flags & (MDOF_X_COMPRESSED)))
 723         {
 724             if (of->tng_low_prec)
 725             {
 726                 real  lambda = -1;
 727                 rvec* box    = nullptr;
 728                 if (mdof_flags & MDOF_BOX_COMPRESSED)
 729                 {
 730                     box = state_local->box;
 731                 }
 732                 if (mdof_flags & MDOF_LAMBDA_COMPRESSED)
 733                 {
 734                     lambda = state_local->lambda[FreeEnergyPerturbationCouplingType::Fep];
 735                 }
 736                 gmx_fwrite_tng(of->tng_low_prec, FALSE, step, t, lambda, box, natoms, nullptr, nullptr, nullptr);
 737             }
 738         }
 739
 740 #if GMX_FAHCORE
 741         /* Write a FAH checkpoint after writing any other data.  We may end up
 742            checkpointing twice but it's fast so it's ok. */
 743         if ((mdof_flags & ~MDOF_CPT))
 744         {
 745             fcCheckpoint();
 746         }
 747 #endif
 748     }
 749 }
 750
 751 void mdoutf_tng_close(gmx_mdoutf_t of)
 752 {
 753     if (of->tng || of->tng_low_prec)
 754     {
 755         wallcycle_start(of->wcycle, ewcTRAJ);
 756         gmx_tng_close(&of->tng);
 757         gmx_tng_close(&of->tng_low_prec);
 758         wallcycle_stop(of->wcycle, ewcTRAJ);
 759     }
 760 }
 761
 762 void done_mdoutf(gmx_mdoutf_t of)
 763 {
 764     if (of->fp_ene != nullptr)
 765     {
 766         done_ener_file(of->fp_ene);
 767     }
 768     if (of->fp_xtc)
 769     {
 770         close_xtc(of->fp_xtc);
 771     }
 772     if (of->fp_trn)
 773     {
 774         gmx_trr_close(of->fp_trn);
 775     }
 776     if (of->fp_dhdl != nullptr)
 777     {
 778         gmx_fio_fclose(of->fp_dhdl);
 779     }
 780     of->outputProvider->finishOutput();
 781     if (of->f_global != nullptr)
 782     {
 783         sfree(of->f_global);
 784     }
 785
 786     gmx_tng_close(&of->tng);
 787     gmx_tng_close(&of->tng_low_prec);
 788
 789     sfree(of);
 790 }
 791
 792 int mdoutf_get_tng_box_output_interval(gmx_mdoutf_t of)
 793 {
 794     if (of->tng)
 795     {
 796         return gmx_tng_get_box_output_interval(of->tng);
 797     }
 798     return 0;
 799 }
 800
 801 int mdoutf_get_tng_lambda_output_interval(gmx_mdoutf_t of)
 802 {
 803     if (of->tng)
 804     {
 805         return gmx_tng_get_lambda_output_interval(of->tng);
 806     }
 807     return 0;
 808 }
 809
 810 int mdoutf_get_tng_compressed_box_output_interval(gmx_mdoutf_t of)
 811 {
 812     if (of->tng_low_prec)
 813     {
 814         return gmx_tng_get_box_output_interval(of->tng_low_prec);
 815     }
 816     return 0;
 817 }
 818
 819 int mdoutf_get_tng_compressed_lambda_output_interval(gmx_mdoutf_t of)
 820 {
 821     if (of->tng_low_prec)
 822     {
 823         return gmx_tng_get_lambda_output_interval(of->tng_low_prec);
 824     }
 825     return 0;
 826 }