src/gromacs/fileio/mdoutf.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2013,2014, 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/fileio/tngio.h"
  40 #include "gromacs/fileio/trajectory_writing.h"
  41 #include "gromacs/fileio/trnio.h"
  42 #include "gromacs/fileio/xtcio.h"
  43 #include "gromacs/fileio/xvgr.h"
  44 #include "gromacs/legacyheaders/checkpoint.h"
  45 #include "gromacs/legacyheaders/copyrite.h"
  46 #include "gromacs/legacyheaders/domdec.h"
  47 #include "gromacs/legacyheaders/mdrun.h"
  48 #include "gromacs/legacyheaders/mvdata.h"
  49 #include "gromacs/legacyheaders/types/commrec.h"
  50 #include "gromacs/math/vec.h"
  51 #include "gromacs/timing/wallcycle.h"
  52 #include "gromacs/utility/fatalerror.h"
  53 #include "gromacs/utility/smalloc.h"
  54
  55 struct gmx_mdoutf {
  56     t_fileio         *fp_trn;
  57     t_fileio         *fp_xtc;
  58     tng_trajectory_t  tng;
  59     tng_trajectory_t  tng_low_prec;
  60     int               x_compression_precision; /* only used by XTC output */
  61     ener_file_t       fp_ene;
  62     const char       *fn_cpt;
  63     gmx_bool          bKeepAndNumCPT;
  64     int               eIntegrator;
  65     gmx_bool          bExpanded;
  66     int               elamstats;
  67     int               simulation_part;
  68     FILE             *fp_dhdl;
  69     FILE             *fp_field;
  70     int               natoms_global;
  71     int               natoms_x_compressed;
  72     gmx_groups_t     *groups; /* for compressed position writing */
  73     gmx_wallcycle_t   wcycle;
  74 };
  75
  76
  77 gmx_mdoutf_t init_mdoutf(FILE *fplog, int nfile, const t_filenm fnm[],
  78                          int mdrun_flags, const t_commrec *cr,
  79                          const t_inputrec *ir, gmx_mtop_t *top_global,
  80                          const output_env_t oenv, gmx_wallcycle_t wcycle)
  81 {
  82     gmx_mdoutf_t  of;
  83     char          filemode[3];
  84     gmx_bool      bAppendFiles, bCiteTng = FALSE;
  85     int           i;
  86
  87     snew(of, 1);
  88
  89     of->fp_trn       = NULL;
  90     of->fp_ene       = NULL;
  91     of->fp_xtc       = NULL;
  92     of->tng          = NULL;
  93     of->tng_low_prec = NULL;
  94     of->fp_dhdl      = NULL;
  95     of->fp_field     = NULL;
  96
  97     of->eIntegrator             = ir->eI;
  98     of->bExpanded               = ir->bExpanded;
  99     of->elamstats               = ir->expandedvals->elamstats;
 100     of->simulation_part         = ir->simulation_part;
 101     of->x_compression_precision = ir->x_compression_precision;
 102     of->wcycle                  = wcycle;
 103
 104     if (MASTER(cr))
 105     {
 106         bAppendFiles = (mdrun_flags & MD_APPENDFILES);
 107
 108         of->bKeepAndNumCPT = (mdrun_flags & MD_KEEPANDNUMCPT);
 109
 110         sprintf(filemode, bAppendFiles ? "a+" : "w+");
 111
 112         if ((EI_DYNAMICS(ir->eI) || EI_ENERGY_MINIMIZATION(ir->eI))
 113 #ifndef GMX_FAHCORE
 114             &&
 115             !(EI_DYNAMICS(ir->eI) &&
 116               ir->nstxout == 0 &&
 117               ir->nstvout == 0 &&
 118               ir->nstfout == 0)
 119 #endif
 120             )
 121         {
 122             const char *filename;
 123             filename = ftp2fn(efTRN, nfile, fnm);
 124             switch (fn2ftp(filename))
 125             {
 126                 case efTRR:
 127                 case efTRN:
 128                     of->fp_trn = open_trn(filename, filemode);
 129                     break;
 130                 case efTNG:
 131                     gmx_tng_open(filename, filemode[0], &of->tng);
 132                     if (filemode[0] == 'w')
 133                     {
 134                         gmx_tng_prepare_md_writing(of->tng, top_global, ir);
 135                     }
 136                     bCiteTng = TRUE;
 137                     break;
 138                 default:
 139                     gmx_incons("Invalid full precision file format");
 140             }
 141         }
 142         if (EI_DYNAMICS(ir->eI) &&
 143             ir->nstxout_compressed > 0)
 144         {
 145             const char *filename;
 146             filename = ftp2fn(efCOMPRESSED, nfile, fnm);
 147             switch (fn2ftp(filename))
 148             {
 149                 case efXTC:
 150                     of->fp_xtc                  = open_xtc(filename, filemode);
 151                     break;
 152                 case efTNG:
 153                     gmx_tng_open(filename, filemode[0], &of->tng_low_prec);
 154                     if (filemode[0] == 'w')
 155                     {
 156                         gmx_tng_prepare_low_prec_writing(of->tng_low_prec, top_global, ir);
 157                     }
 158                     bCiteTng = TRUE;
 159                     break;
 160                 default:
 161                     gmx_incons("Invalid reduced precision file format");
 162             }
 163         }
 164         if (EI_DYNAMICS(ir->eI) || EI_ENERGY_MINIMIZATION(ir->eI))
 165         {
 166             of->fp_ene = open_enx(ftp2fn(efEDR, nfile, fnm), filemode);
 167         }
 168         of->fn_cpt = opt2fn("-cpo", nfile, fnm);
 169
 170         if ((ir->efep != efepNO || ir->bSimTemp) && ir->fepvals->nstdhdl > 0 &&
 171             (ir->fepvals->separate_dhdl_file == esepdhdlfileYES ) &&
 172             EI_DYNAMICS(ir->eI))
 173         {
 174             if (bAppendFiles)
 175             {
 176                 of->fp_dhdl = gmx_fio_fopen(opt2fn("-dhdl", nfile, fnm), filemode);
 177             }
 178             else
 179             {
 180                 of->fp_dhdl = open_dhdl(opt2fn("-dhdl", nfile, fnm), ir, oenv);
 181             }
 182         }
 183
 184         if (opt2bSet("-field", nfile, fnm) &&
 185             (ir->ex[XX].n || ir->ex[YY].n || ir->ex[ZZ].n))
 186         {
 187             if (bAppendFiles)
 188             {
 189                 of->fp_dhdl = gmx_fio_fopen(opt2fn("-field", nfile, fnm),
 190                                             filemode);
 191             }
 192             else
 193             {
 194                 of->fp_field = xvgropen(opt2fn("-field", nfile, fnm),
 195                                         "Applied electric field", "Time (ps)",
 196                                         "E (V/nm)", oenv);
 197             }
 198         }
 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 (ggrpnr(of->groups, egcCompressedX, i) == 0)
 210             {
 211                 of->natoms_x_compressed++;
 212             }
 213         }
 214     }
 215
 216     if (bCiteTng)
 217     {
 218         please_cite(fplog, "Lundborg2014");
 219     }
 220
 221     return of;
 222 }
 223
 224 FILE *mdoutf_get_fp_field(gmx_mdoutf_t of)
 225 {
 226     return of->fp_field;
 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, t_commrec *cr,
 245                                       gmx_mdoutf_t of,
 246                                       int mdof_flags,
 247                                       gmx_mtop_t *top_global,
 248                                       gmx_int64_t step, double t,
 249                                       t_state *state_local, t_state *state_global,
 250                                       rvec *f_local, rvec *f_global)
 251 {
 252     rvec *local_v;
 253     rvec *global_v;
 254
 255     /* MRS -- defining these variables is to manage the difference
 256      * between half step and full step velocities, but there must be a better way . . . */
 257
 258     local_v  = state_local->v;
 259     global_v = state_global->v;
 260
 261     if (DOMAINDECOMP(cr))
 262     {
 263         if (mdof_flags & MDOF_CPT)
 264         {
 265             dd_collect_state(cr->dd, state_local, state_global);
 266         }
 267         else
 268         {
 269             if (mdof_flags & (MDOF_X | MDOF_X_COMPRESSED))
 270             {
 271                 dd_collect_vec(cr->dd, state_local, state_local->x,
 272                                state_global->x);
 273             }
 274             if (mdof_flags & MDOF_V)
 275             {
 276                 dd_collect_vec(cr->dd, state_local, local_v,
 277                                global_v);
 278             }
 279         }
 280         if (mdof_flags & MDOF_F)
 281         {
 282             dd_collect_vec(cr->dd, state_local, f_local, f_global);
 283         }
 284     }
 285     else
 286     {
 287         if (mdof_flags & MDOF_CPT)
 288         {
 289             /* All pointers in state_local are equal to state_global,
 290              * but we need to copy the non-pointer entries.
 291              */
 292             state_global->lambda = state_local->lambda;
 293             state_global->veta   = state_local->veta;
 294             state_global->vol0   = state_local->vol0;
 295             copy_mat(state_local->box, state_global->box);
 296             copy_mat(state_local->boxv, state_global->boxv);
 297             copy_mat(state_local->svir_prev, state_global->svir_prev);
 298             copy_mat(state_local->fvir_prev, state_global->fvir_prev);
 299             copy_mat(state_local->pres_prev, state_global->pres_prev);
 300         }
 301     }
 302
 303     if (MASTER(cr))
 304     {
 305         if (mdof_flags & MDOF_CPT)
 306         {
 307             fflush_tng(of->tng);
 308             fflush_tng(of->tng_low_prec);
 309             write_checkpoint(of->fn_cpt, of->bKeepAndNumCPT,
 310                              fplog, cr, of->eIntegrator, of->simulation_part,
 311                              of->bExpanded, of->elamstats, step, t, state_global);
 312         }
 313
 314         if (mdof_flags & (MDOF_X | MDOF_V | MDOF_F))
 315         {
 316             if (of->fp_trn)
 317             {
 318                 fwrite_trn(of->fp_trn, step, t, state_local->lambda[efptFEP],
 319                            state_local->box, top_global->natoms,
 320                            (mdof_flags & MDOF_X) ? state_global->x : NULL,
 321                            (mdof_flags & MDOF_V) ? global_v : NULL,
 322                            (mdof_flags & MDOF_F) ? f_global : NULL);
 323                 if (gmx_fio_flush(of->fp_trn) != 0)
 324                 {
 325                     gmx_file("Cannot write trajectory; maybe you are out of disk space?");
 326                 }
 327             }
 328
 329             gmx_fwrite_tng(of->tng, FALSE, step, t, state_local->lambda[efptFEP],
 330                            (const rvec *) state_local->box,
 331                            top_global->natoms,
 332                            (mdof_flags & MDOF_X) ? (const rvec *) state_global->x : NULL,
 333                            (mdof_flags & MDOF_V) ? (const rvec *) global_v : NULL,
 334                            (mdof_flags & MDOF_F) ? (const rvec *) f_global : NULL);
 335         }
 336         if (mdof_flags & MDOF_X_COMPRESSED)
 337         {
 338             rvec *xxtc = NULL;
 339
 340             if (of->natoms_x_compressed == of->natoms_global)
 341             {
 342                 /* We are writing the positions of all of the atoms to
 343                    the compressed output */
 344                 xxtc = state_global->x;
 345             }
 346             else
 347             {
 348                 /* We are writing the positions of only a subset of
 349                    the atoms to the compressed output, so we have to
 350                    make a copy of the subset of coordinates. */
 351                 int i, j;
 352
 353                 snew(xxtc, of->natoms_x_compressed);
 354                 for (i = 0, j = 0; (i < of->natoms_global); i++)
 355                 {
 356                     if (ggrpnr(of->groups, egcCompressedX, i) == 0)
 357                     {
 358                         copy_rvec(state_global->x[i], xxtc[j++]);
 359                     }
 360                 }
 361             }
 362             if (write_xtc(of->fp_xtc, of->natoms_x_compressed, step, t,
 363                           state_local->box, xxtc, of->x_compression_precision) == 0)
 364             {
 365                 gmx_fatal(FARGS, "XTC error - maybe you are out of disk space?");
 366             }
 367             gmx_fwrite_tng(of->tng_low_prec,
 368                            TRUE,
 369                            step,
 370                            t,
 371                            state_local->lambda[efptFEP],
 372                            (const rvec *) state_local->box,
 373                            of->natoms_x_compressed,
 374                            (const rvec *) xxtc,
 375                            NULL,
 376                            NULL);
 377             if (of->natoms_x_compressed != of->natoms_global)
 378             {
 379                 sfree(xxtc);
 380             }
 381         }
 382     }
 383 }
 384
 385 void mdoutf_tng_close(gmx_mdoutf_t of)
 386 {
 387     if (of->tng || of->tng_low_prec)
 388     {
 389         wallcycle_start(of->wcycle, ewcTRAJ);
 390         gmx_tng_close(&of->tng);
 391         gmx_tng_close(&of->tng_low_prec);
 392         wallcycle_stop(of->wcycle, ewcTRAJ);
 393     }
 394 }
 395
 396 void done_mdoutf(gmx_mdoutf_t of)
 397 {
 398     if (of->fp_ene != NULL)
 399     {
 400         close_enx(of->fp_ene);
 401     }
 402     if (of->fp_xtc)
 403     {
 404         close_xtc(of->fp_xtc);
 405     }
 406     if (of->fp_trn)
 407     {
 408         close_trn(of->fp_trn);
 409     }
 410     if (of->fp_dhdl != NULL)
 411     {
 412         gmx_fio_fclose(of->fp_dhdl);
 413     }
 414     if (of->fp_field != NULL)
 415     {
 416         gmx_fio_fclose(of->fp_field);
 417     }
 418
 419     gmx_tng_close(&of->tng);
 420     gmx_tng_close(&of->tng_low_prec);
 421
 422     sfree(of);
 423 }