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