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