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