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