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