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