src/gromacs/fileio/mdoutf.c

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