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