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