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