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