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