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/smalloc.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
  48 struct gmx_mdoutf {
  49     t_fileio         *fp_trn;
  50     t_fileio         *fp_xtc;
  51     tng_trajectory_t  tng;
  52     tng_trajectory_t  tng_low_prec;
  53     int               x_compression_precision; /* only used by XTC output */
  54     ener_file_t       fp_ene;
  55     const char       *fn_cpt;
  56     gmx_bool          bKeepAndNumCPT;
  57     int               eIntegrator;
  58     gmx_bool          bExpanded;
  59     int               elamstats;
  60     int               simulation_part;
  61     FILE             *fp_dhdl;
  62     FILE             *fp_field;
  63     int               natoms_global;
  64     int               natoms_x_compressed;
  65     gmx_groups_t     *groups; /* for compressed position writing */
  66 };
  67
  68
  69 gmx_mdoutf_t init_mdoutf(int nfile, const t_filenm fnm[], int mdrun_flags,
  70                          const t_commrec *cr, const t_inputrec *ir,
  71                          gmx_mtop_t *top_global,
  72                          const output_env_t oenv)
  73 {
  74     gmx_mdoutf_t  of;
  75     char          filemode[3];
  76     gmx_bool      bAppendFiles;
  77     int           i;
  78
  79     snew(of, 1);
  80
  81     of->fp_trn       = NULL;
  82     of->fp_ene       = NULL;
  83     of->fp_xtc       = NULL;
  84     of->tng          = NULL;
  85     of->tng_low_prec = NULL;
  86     of->fp_dhdl      = NULL;
  87     of->fp_field     = NULL;
  88
  89     of->eIntegrator             = ir->eI;
  90     of->bExpanded               = ir->bExpanded;
  91     of->elamstats               = ir->expandedvals->elamstats;
  92     of->simulation_part         = ir->simulation_part;
  93     of->x_compression_precision = ir->x_compression_precision;
  94
  95     if (MASTER(cr))
  96     {
  97         bAppendFiles = (mdrun_flags & MD_APPENDFILES);
  98
  99         of->bKeepAndNumCPT = (mdrun_flags & MD_KEEPANDNUMCPT);
 100
 101         sprintf(filemode, bAppendFiles ? "a+" : "w+");
 102
 103         if ((EI_DYNAMICS(ir->eI) || EI_ENERGY_MINIMIZATION(ir->eI))
 104 #ifndef GMX_FAHCORE
 105             &&
 106             !(EI_DYNAMICS(ir->eI) &&
 107               ir->nstxout == 0 &&
 108               ir->nstvout == 0 &&
 109               ir->nstfout == 0)
 110 #endif
 111             )
 112         {
 113             const char *filename;
 114             filename = ftp2fn(efTRN, nfile, fnm);
 115             switch (fn2ftp(filename))
 116             {
 117                 case efTRR:
 118                 case efTRN:
 119                     of->fp_trn = open_trn(filename, filemode);
 120                     break;
 121                 case efTNG:
 122                     gmx_tng_open(filename, filemode[0], &of->tng);
 123                     if (filemode[0] == 'w')
 124                     {
 125                         gmx_tng_prepare_md_writing(of->tng, top_global, ir);
 126                     }
 127                     break;
 128                 default:
 129                     gmx_incons("Invalid full precision file format");
 130             }
 131         }
 132         if (EI_DYNAMICS(ir->eI) &&
 133             ir->nstxout_compressed > 0)
 134         {
 135             const char *filename;
 136             filename = ftp2fn(efCOMPRESSED, nfile, fnm);
 137             switch (fn2ftp(filename))
 138             {
 139                 case efXTC:
 140                     of->fp_xtc                  = open_xtc(filename, filemode);
 141                     break;
 142                 case efTNG:
 143                     gmx_tng_open(filename, filemode[0], &of->tng_low_prec);
 144                     if (filemode[0] == 'w')
 145                     {
 146                         gmx_tng_prepare_low_prec_writing(of->tng_low_prec, top_global, ir);
 147                     }
 148                     break;
 149                 default:
 150                     gmx_incons("Invalid reduced precision file format");
 151             }
 152         }
 153         if (EI_DYNAMICS(ir->eI) || EI_ENERGY_MINIMIZATION(ir->eI))
 154         {
 155             of->fp_ene = open_enx(ftp2fn(efEDR, nfile, fnm), filemode);
 156         }
 157         of->fn_cpt = opt2fn("-cpo", nfile, fnm);
 158
 159         if ((ir->efep != efepNO || ir->bSimTemp) && ir->fepvals->nstdhdl > 0 &&
 160             (ir->fepvals->separate_dhdl_file == esepdhdlfileYES ) &&
 161             EI_DYNAMICS(ir->eI))
 162         {
 163             if (bAppendFiles)
 164             {
 165                 of->fp_dhdl = gmx_fio_fopen(opt2fn("-dhdl", nfile, fnm), filemode);
 166             }
 167             else
 168             {
 169                 of->fp_dhdl = open_dhdl(opt2fn("-dhdl", nfile, fnm), ir, oenv);
 170             }
 171         }
 172
 173         if (opt2bSet("-field", nfile, fnm) &&
 174             (ir->ex[XX].n || ir->ex[YY].n || ir->ex[ZZ].n))
 175         {
 176             if (bAppendFiles)
 177             {
 178                 of->fp_dhdl = gmx_fio_fopen(opt2fn("-field", nfile, fnm),
 179                                             filemode);
 180             }
 181             else
 182             {
 183                 of->fp_field = xvgropen(opt2fn("-field", nfile, fnm),
 184                                         "Applied electric field", "Time (ps)",
 185                                         "E (V/nm)", oenv);
 186             }
 187         }
 188
 189         /* Set up atom counts so they can be passed to actual
 190            trajectory-writing routines later. Also, XTC writing needs
 191            to know what (and how many) atoms might be in the XTC
 192            groups, and how to look up later which ones they are. */
 193         of->natoms_global       = top_global->natoms;
 194         of->groups              = &top_global->groups;
 195         of->natoms_x_compressed = 0;
 196         for (i = 0; (i < top_global->natoms); i++)
 197         {
 198             if (ggrpnr(of->groups, egcCompressedX, i) == 0)
 199             {
 200                 of->natoms_x_compressed++;
 201             }
 202         }
 203     }
 204
 205     return of;
 206 }
 207
 208 FILE *mdoutf_get_fp_field(gmx_mdoutf_t of)
 209 {
 210     return of->fp_field;
 211 }
 212
 213 ener_file_t mdoutf_get_fp_ene(gmx_mdoutf_t of)
 214 {
 215     return of->fp_ene;
 216 }
 217
 218 FILE *mdoutf_get_fp_dhdl(gmx_mdoutf_t of)
 219 {
 220     return of->fp_dhdl;
 221 }
 222
 223 static void moveit(t_commrec *cr, rvec xx[])
 224 {
 225     if (!xx)
 226     {
 227         return;
 228     }
 229
 230     move_rvecs(cr, FALSE, FALSE, xx, NULL, (cr->nnodes-cr->npmenodes)-1, NULL);
 231 }
 232
 233 void mdoutf_write_to_trajectory_files(FILE *fplog, t_commrec *cr,
 234                                       gmx_mdoutf_t of,
 235                                       int mdof_flags,
 236                                       gmx_mtop_t *top_global,
 237                                       gmx_int64_t step, double t,
 238                                       t_state *state_local, t_state *state_global,
 239                                       rvec *f_local, rvec *f_global)
 240 {
 241     rvec *local_v;
 242     rvec *global_v;
 243
 244 #define MX(xvf) moveit(cr, xvf)
 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         if (cr->nnodes > 1)
 293         {
 294             /* Particle decomposition, collect the data on the master node */
 295             if (mdof_flags & MDOF_CPT)
 296             {
 297                 if (state_local->flags & (1<<estX))
 298                 {
 299                     MX(state_global->x);
 300                 }
 301                 if (state_local->flags & (1<<estV))
 302                 {
 303                     MX(state_global->v);
 304                 }
 305                 if (state_local->flags & (1<<estSDX))
 306                 {
 307                     MX(state_global->sd_X);
 308                 }
 309                 if (state_global->nrngi > 1)
 310                 {
 311                     if (state_local->flags & (1<<estLD_RNG))
 312                     {
 313 #ifdef GMX_MPI
 314                         MPI_Gather(state_local->ld_rng,
 315                                    state_local->nrng*sizeof(state_local->ld_rng[0]), MPI_BYTE,
 316                                    state_global->ld_rng,
 317                                    state_local->nrng*sizeof(state_local->ld_rng[0]), MPI_BYTE,
 318                                    MASTERRANK(cr), cr->mpi_comm_mygroup);
 319 #endif
 320                     }
 321                     if (state_local->flags & (1<<estLD_RNGI))
 322                     {
 323 #ifdef GMX_MPI
 324                         MPI_Gather(state_local->ld_rngi,
 325                                    sizeof(state_local->ld_rngi[0]), MPI_BYTE,
 326                                    state_global->ld_rngi,
 327                                    sizeof(state_local->ld_rngi[0]), MPI_BYTE,
 328                                    MASTERRANK(cr), cr->mpi_comm_mygroup);
 329 #endif
 330                     }
 331                 }
 332             }
 333             else
 334             {
 335                 if (mdof_flags & (MDOF_X | MDOF_X_COMPRESSED))
 336                 {
 337                     MX(state_global->x);
 338                 }
 339                 if (mdof_flags & MDOF_V)
 340                 {
 341                     MX(global_v);
 342                 }
 343             }
 344             if (mdof_flags & MDOF_F)
 345             {
 346                 MX(f_global);
 347             }
 348         }
 349     }
 350
 351     if (MASTER(cr))
 352     {
 353         if (mdof_flags & MDOF_CPT)
 354         {
 355             fflush_tng(of->tng);
 356             fflush_tng(of->tng_low_prec);
 357             write_checkpoint(of->fn_cpt, of->bKeepAndNumCPT,
 358                              fplog, cr, of->eIntegrator, of->simulation_part,
 359                              of->bExpanded, of->elamstats, step, t, state_global);
 360         }
 361
 362         if (mdof_flags & (MDOF_X | MDOF_V | MDOF_F))
 363         {
 364             if (of->fp_trn)
 365             {
 366                 fwrite_trn(of->fp_trn, step, t, state_local->lambda[efptFEP],
 367                            state_local->box, top_global->natoms,
 368                            (mdof_flags & MDOF_X) ? state_global->x : NULL,
 369                            (mdof_flags & MDOF_V) ? global_v : NULL,
 370                            (mdof_flags & MDOF_F) ? f_global : NULL);
 371                 if (gmx_fio_flush(of->fp_trn) != 0)
 372                 {
 373                     gmx_file("Cannot write trajectory; maybe you are out of disk space?");
 374                 }
 375             }
 376
 377             gmx_fwrite_tng(of->tng, FALSE, step, t, state_local->lambda[efptFEP],
 378                            (const rvec *) state_local->box,
 379                            top_global->natoms,
 380                            (mdof_flags & MDOF_X) ? (const rvec *) state_global->x : NULL,
 381                            (mdof_flags & MDOF_V) ? (const rvec *) global_v : NULL,
 382                            (mdof_flags & MDOF_F) ? (const rvec *) f_global : NULL);
 383         }
 384         if (mdof_flags & MDOF_X_COMPRESSED)
 385         {
 386             rvec *xxtc = NULL;
 387
 388             if (of->natoms_x_compressed == of->natoms_global)
 389             {
 390                 /* We are writing the positions of all of the atoms to
 391                    the compressed output */
 392                 xxtc = state_global->x;
 393             }
 394             else
 395             {
 396                 /* We are writing the positions of only a subset of
 397                    the atoms to the compressed output, so we have to
 398                    make a copy of the subset of coordinates. */
 399                 int i, j;
 400
 401                 snew(xxtc, of->natoms_x_compressed);
 402                 for (i = 0, j = 0; (i < of->natoms_x_compressed); i++)
 403                 {
 404                     if (ggrpnr(of->groups, egcCompressedX, i) == 0)
 405                     {
 406                         copy_rvec(state_global->x[i], xxtc[j++]);
 407                     }
 408                 }
 409             }
 410             if (write_xtc(of->fp_xtc, of->natoms_x_compressed, step, t,
 411                           state_local->box, xxtc, of->x_compression_precision) == 0)
 412             {
 413                 gmx_fatal(FARGS, "XTC error - maybe you are out of disk space?");
 414             }
 415             gmx_fwrite_tng(of->tng_low_prec,
 416                            TRUE,
 417                            step,
 418                            t,
 419                            state_local->lambda[efptFEP],
 420                            (const rvec *) state_local->box,
 421                            of->natoms_x_compressed,
 422                            (const rvec *) xxtc,
 423                            NULL,
 424                            NULL);
 425             if (of->natoms_x_compressed != of->natoms_global)
 426             {
 427                 sfree(xxtc);
 428             }
 429         }
 430     }
 431 }
 432
 433 void done_mdoutf(gmx_mdoutf_t of)
 434 {
 435     if (of->fp_ene != NULL)
 436     {
 437         close_enx(of->fp_ene);
 438     }
 439     if (of->fp_xtc)
 440     {
 441         close_xtc(of->fp_xtc);
 442     }
 443     if (of->fp_trn)
 444     {
 445         close_trn(of->fp_trn);
 446     }
 447     if (of->fp_dhdl != NULL)
 448     {
 449         gmx_fio_fclose(of->fp_dhdl);
 450     }
 451     if (of->fp_field != NULL)
 452     {
 453         gmx_fio_fclose(of->fp_field);
 454     }
 455     gmx_tng_close(&of->tng);
 456     gmx_tng_close(&of->tng_low_prec);
 457
 458     sfree(of);
 459 }