src/mdlib/stat.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <string.h>
  43 #include <stdio.h>
  44 #include "typedefs.h"
  45 #include "sysstuff.h"
  46 #include "gmx_fatal.h"
  47 #include "network.h"
  48 #include "txtdump.h"
  49 #include "names.h"
  50 #include "physics.h"
  51 #include "vec.h"
  52 #include "maths.h"
  53 #include "mvdata.h"
  54 #include "main.h"
  55 #include "force.h"
  56 #include "vcm.h"
  57 #include "smalloc.h"
  58 #include "futil.h"
  59 #include "network.h"
  60 #include "rbin.h"
  61 #include "tgroup.h"
  62 #include "xtcio.h"
  63 #include "gmxfio.h"
  64 #include "trnio.h"
  65 #include "statutil.h"
  66 #include "domdec.h"
  67 #include "partdec.h"
  68 #include "constr.h"
  69 #include "checkpoint.h"
  70 #include "xvgr.h"
  71 #include "md_support.h"
  72 #include "mdrun.h"
  73 #include "sim_util.h"
  74
  75 typedef struct gmx_global_stat
  76 {
  77     t_bin *rb;
  78     int   *itc0;
  79     int   *itc1;
  80 } t_gmx_global_stat;
  81
  82 gmx_global_stat_t global_stat_init(t_inputrec *ir)
  83 {
  84     gmx_global_stat_t gs;
  85
  86     snew(gs,1);
  87
  88     gs->rb = mk_bin();
  89     snew(gs->itc0,ir->opts.ngtc);
  90     snew(gs->itc1,ir->opts.ngtc);
  91
  92     return gs;
  93 }
  94
  95 void global_stat_destroy(gmx_global_stat_t gs)
  96 {
  97     destroy_bin(gs->rb);
  98     sfree(gs->itc0);
  99     sfree(gs->itc1);
 100     sfree(gs);
 101 }
 102
 103 static int filter_enerdterm(real *afrom, gmx_bool bToBuffer, real *ato,
 104                             gmx_bool bTemp, gmx_bool bPres, gmx_bool bEner) {
 105     int i,to,from;
 106
 107     from = 0;
 108     to   = 0;
 109     for (i=0;i<F_NRE;i++)
 110     {
 111         if (bToBuffer)
 112         {
 113             from = i;
 114         }
 115         else
 116         {
 117             to = i;
 118         }
 119         switch (i) {
 120         case F_EKIN:
 121         case F_TEMP:
 122         case F_DKDL:
 123             if (bTemp)
 124             {
 125                 ato[to++] = afrom[from++];
 126             }
 127             break;
 128         case F_PRES:
 129         case F_PDISPCORR:
 130         case F_VTEMP:
 131             if (bPres)
 132             {
 133                 ato[to++] = afrom[from++];
 134             }
 135             break;
 136         default:
 137             if (bEner)
 138             {
 139                 ato[to++] = afrom[from++];
 140             }
 141             break;
 142         }
 143     }
 144
 145     return to;
 146 }
 147
 148 void global_stat(FILE *fplog,gmx_global_stat_t gs,
 149                  t_commrec *cr,gmx_enerdata_t *enerd,
 150                  tensor fvir,tensor svir,rvec mu_tot,
 151                  t_inputrec *inputrec,
 152                  gmx_ekindata_t *ekind,gmx_constr_t constr,
 153                  t_vcm *vcm,
 154                  int nsig,real *sig,
 155                  gmx_mtop_t *top_global, t_state *state_local,
 156                  gmx_bool bSumEkinhOld, int flags)
 157 /* instead of current system, gmx_booleans for summing virial, kinetic energy, and other terms */
 158 {
 159   t_bin  *rb;
 160   int    *itc0,*itc1;
 161   int    ie=0,ifv=0,isv=0,irmsd=0,imu=0;
 162   int    idedl=0,idvdll=0,idvdlnl=0,iepl=0,icm=0,imass=0,ica=0,inb=0;
 163   int    isig=-1;
 164   int    icj=-1,ici=-1,icx=-1;
 165   int    inn[egNR];
 166   real   copyenerd[F_NRE];
 167   int    nener,j;
 168   real   *rmsd_data=NULL;
 169   double nb;
 170   gmx_bool   bVV,bTemp,bEner,bPres,bConstrVir,bEkinAveVel,bFirstIterate,bReadEkin;
 171
 172   bVV           = EI_VV(inputrec->eI);
 173   bTemp         = flags & CGLO_TEMPERATURE;
 174   bEner         = flags & CGLO_ENERGY;
 175   bPres         = (flags & CGLO_PRESSURE);
 176   bConstrVir    = (flags & CGLO_CONSTRAINT);
 177   bFirstIterate = (flags & CGLO_FIRSTITERATE);
 178   bEkinAveVel   = (inputrec->eI==eiVV || (inputrec->eI==eiVVAK && bPres));
 179   bReadEkin     = (flags & CGLO_READEKIN);
 180
 181   rb   = gs->rb;
 182   itc0 = gs->itc0;
 183   itc1 = gs->itc1;
 184
 185
 186   reset_bin(rb);
 187   /* This routine copies all the data to be summed to one big buffer
 188    * using the t_bin struct.
 189    */
 190
 191   /* First, we neeed to identify which enerd->term should be
 192      communicated.  Temperature and pressure terms should only be
 193      communicated and summed when they need to be, to avoid repeating
 194      the sums and overcounting. */
 195
 196   nener = filter_enerdterm(enerd->term,TRUE,copyenerd,bTemp,bPres,bEner);
 197
 198   /* First, the data that needs to be communicated with velocity verlet every time
 199      This is just the constraint virial.*/
 200   if (bConstrVir) {
 201       isv = add_binr(rb,DIM*DIM,svir[0]);
 202       where();
 203   }
 204
 205 /* We need the force virial and the kinetic energy for the first time through with velocity verlet */
 206   if (bTemp || !bVV)
 207   {
 208       if (ekind)
 209       {
 210           for(j=0; (j<inputrec->opts.ngtc); j++)
 211           {
 212               if (bSumEkinhOld)
 213               {
 214                   itc0[j]=add_binr(rb,DIM*DIM,ekind->tcstat[j].ekinh_old[0]);
 215               }
 216               if (bEkinAveVel && !bReadEkin)
 217               {
 218                   itc1[j]=add_binr(rb,DIM*DIM,ekind->tcstat[j].ekinf[0]);
 219               }
 220               else if (!bReadEkin)
 221               {
 222                   itc1[j]=add_binr(rb,DIM*DIM,ekind->tcstat[j].ekinh[0]);
 223               }
 224           }
 225           /* these probably need to be put into one of these categories */
 226           where();
 227           idedl = add_binr(rb,1,&(ekind->dekindl));
 228           where();
 229           ica   = add_binr(rb,1,&(ekind->cosacc.mvcos));
 230           where();
 231       }
 232   }
 233   where();
 234
 235   if ((bPres || !bVV) && bFirstIterate)
 236   {
 237       ifv = add_binr(rb,DIM*DIM,fvir[0]);
 238   }
 239
 240
 241   if (bEner)
 242   {
 243       where();
 244       if (bFirstIterate)
 245       {
 246           ie  = add_binr(rb,nener,copyenerd);
 247       }
 248       where();
 249       if (constr)
 250       {
 251           rmsd_data = constr_rmsd_data(constr);
 252           if (rmsd_data)
 253           {
 254               irmsd = add_binr(rb,inputrec->eI==eiSD2 ? 3 : 2,rmsd_data);
 255           }
 256       }
 257       if (!NEED_MUTOT(*inputrec))
 258       {
 259           imu = add_binr(rb,DIM,mu_tot);
 260           where();
 261       }
 262
 263       if (bFirstIterate)
 264       {
 265           for(j=0; (j<egNR); j++)
 266           {
 267               inn[j]=add_binr(rb,enerd->grpp.nener,enerd->grpp.ener[j]);
 268           }
 269           where();
 270           if (inputrec->efep != efepNO)
 271           {
 272               idvdll  = add_bind(rb,efptNR,enerd->dvdl_lin);
 273               idvdlnl = add_bind(rb,efptNR,enerd->dvdl_nonlin);
 274               if (enerd->n_lambda > 0)
 275               {
 276                   iepl = add_bind(rb,enerd->n_lambda,enerd->enerpart_lambda);
 277               }
 278           }
 279       }
 280   }
 281
 282   if (vcm)
 283   {
 284       icm   = add_binr(rb,DIM*vcm->nr,vcm->group_p[0]);
 285       where();
 286       imass = add_binr(rb,vcm->nr,vcm->group_mass);
 287       where();
 288       if (vcm->mode == ecmANGULAR)
 289       {
 290           icj   = add_binr(rb,DIM*vcm->nr,vcm->group_j[0]);
 291           where();
 292           icx   = add_binr(rb,DIM*vcm->nr,vcm->group_x[0]);
 293           where();
 294           ici   = add_binr(rb,DIM*DIM*vcm->nr,vcm->group_i[0][0]);
 295           where();
 296       }
 297   }
 298
 299   if (DOMAINDECOMP(cr))
 300   {
 301       nb = cr->dd->nbonded_local;
 302       inb = add_bind(rb,1,&nb);
 303       }
 304   where();
 305   if (nsig > 0)
 306   {
 307       isig = add_binr(rb,nsig,sig);
 308   }
 309
 310   /* Global sum it all */
 311   if (debug)
 312   {
 313       fprintf(debug,"Summing %d energies\n",rb->maxreal);
 314   }
 315   sum_bin(rb,cr);
 316   where();
 317
 318   /* Extract all the data locally */
 319
 320   if (bConstrVir)
 321   {
 322       extract_binr(rb,isv ,DIM*DIM,svir[0]);
 323   }
 324
 325   /* We need the force virial and the kinetic energy for the first time through with velocity verlet */
 326   if (bTemp || !bVV)
 327   {
 328       if (ekind)
 329       {
 330           for(j=0; (j<inputrec->opts.ngtc); j++)
 331           {
 332               if (bSumEkinhOld)
 333               {
 334                   extract_binr(rb,itc0[j],DIM*DIM,ekind->tcstat[j].ekinh_old[0]);
 335               }
 336               if (bEkinAveVel && !bReadEkin) {
 337                   extract_binr(rb,itc1[j],DIM*DIM,ekind->tcstat[j].ekinf[0]);
 338               }
 339               else if (!bReadEkin)
 340               {
 341                   extract_binr(rb,itc1[j],DIM*DIM,ekind->tcstat[j].ekinh[0]);
 342               }
 343           }
 344           extract_binr(rb,idedl,1,&(ekind->dekindl));
 345           extract_binr(rb,ica,1,&(ekind->cosacc.mvcos));
 346           where();
 347       }
 348   }
 349   if ((bPres || !bVV) && bFirstIterate)
 350   {
 351       extract_binr(rb,ifv ,DIM*DIM,fvir[0]);
 352   }
 353
 354   if (bEner)
 355   {
 356       if (bFirstIterate)
 357       {
 358           extract_binr(rb,ie,nener,copyenerd);
 359           if (rmsd_data)
 360           {
 361               extract_binr(rb,irmsd,inputrec->eI==eiSD2 ? 3 : 2,rmsd_data);
 362           }
 363           if (!NEED_MUTOT(*inputrec))
 364           {
 365               extract_binr(rb,imu,DIM,mu_tot);
 366           }
 367
 368           for(j=0; (j<egNR); j++)
 369           {
 370               extract_binr(rb,inn[j],enerd->grpp.nener,enerd->grpp.ener[j]);
 371           }
 372           if (inputrec->efep != efepNO)
 373           {
 374               extract_bind(rb,idvdll ,efptNR,enerd->dvdl_lin);
 375               extract_bind(rb,idvdlnl,efptNR,enerd->dvdl_nonlin);
 376               if (enerd->n_lambda > 0)
 377               {
 378                   extract_bind(rb,iepl,enerd->n_lambda,enerd->enerpart_lambda);
 379               }
 380           }
 381           if (DOMAINDECOMP(cr))
 382           {
 383               extract_bind(rb,inb,1,&nb);
 384               if ((int)(nb + 0.5) != cr->dd->nbonded_global)
 385               {
 386                   dd_print_missing_interactions(fplog,cr,(int)(nb + 0.5),top_global,state_local);
 387               }
 388           }
 389           where();
 390
 391           filter_enerdterm(copyenerd,FALSE,enerd->term,bTemp,bPres,bEner);
 392       }
 393   }
 394
 395   if (vcm)
 396   {
 397       extract_binr(rb,icm,DIM*vcm->nr,vcm->group_p[0]);
 398       where();
 399       extract_binr(rb,imass,vcm->nr,vcm->group_mass);
 400       where();
 401       if (vcm->mode == ecmANGULAR)
 402       {
 403           extract_binr(rb,icj,DIM*vcm->nr,vcm->group_j[0]);
 404           where();
 405           extract_binr(rb,icx,DIM*vcm->nr,vcm->group_x[0]);
 406           where();
 407           extract_binr(rb,ici,DIM*DIM*vcm->nr,vcm->group_i[0][0]);
 408           where();
 409       }
 410   }
 411
 412   if (nsig > 0)
 413   {
 414       extract_binr(rb,isig,nsig,sig);
 415   }
 416   where();
 417 }
 418
 419 int do_per_step(gmx_large_int_t step,gmx_large_int_t nstep)
 420 {
 421   if (nstep != 0)
 422     return ((step % nstep)==0);
 423   else
 424     return 0;
 425 }
 426
 427 static void moveit(t_commrec *cr,
 428                    int left,int right,const char *s,rvec xx[])
 429 {
 430   if (!xx)
 431     return;
 432
 433   move_rvecs(cr,FALSE,FALSE,left,right,
 434              xx,NULL,(cr->nnodes-cr->npmenodes)-1,NULL);
 435 }
 436
 437 gmx_mdoutf_t *init_mdoutf(int nfile,const t_filenm fnm[],int mdrun_flags,
 438                           const t_commrec *cr,const t_inputrec *ir,
 439                           const output_env_t oenv)
 440 {
 441     gmx_mdoutf_t *of;
 442     char filemode[3];
 443     gmx_bool bAppendFiles;
 444
 445     snew(of,1);
 446
 447     of->fp_trn   = NULL;
 448     of->fp_ene   = NULL;
 449     of->fp_xtc   = NULL;
 450     of->fp_dhdl  = NULL;
 451     of->fp_field = NULL;
 452
 453     of->eIntegrator     = ir->eI;
 454     of->bExpanded       = ir->bExpanded;
 455     of->elamstats       = ir->expandedvals->elamstats;
 456     of->simulation_part = ir->simulation_part;
 457
 458     if (MASTER(cr))
 459     {
 460         bAppendFiles = (mdrun_flags & MD_APPENDFILES);
 461
 462         of->bKeepAndNumCPT = (mdrun_flags & MD_KEEPANDNUMCPT);
 463
 464         sprintf(filemode, bAppendFiles ? "a+" : "w+");
 465
 466         if ((EI_DYNAMICS(ir->eI) || EI_ENERGY_MINIMIZATION(ir->eI))
 467 #ifndef GMX_FAHCORE
 468             &&
 469             !(EI_DYNAMICS(ir->eI) &&
 470               ir->nstxout == 0 &&
 471               ir->nstvout == 0 &&
 472               ir->nstfout == 0)
 473 #endif
 474             )
 475         {
 476             of->fp_trn = open_trn(ftp2fn(efTRN,nfile,fnm), filemode);
 477         }
 478         if (EI_DYNAMICS(ir->eI) &&
 479             ir->nstxtcout > 0)
 480         {
 481             of->fp_xtc = open_xtc(ftp2fn(efXTC,nfile,fnm), filemode);
 482             of->xtc_prec = ir->xtcprec;
 483         }
 484         if (EI_DYNAMICS(ir->eI) || EI_ENERGY_MINIMIZATION(ir->eI))
 485         {
 486             of->fp_ene = open_enx(ftp2fn(efEDR,nfile,fnm), filemode);
 487         }
 488         of->fn_cpt = opt2fn("-cpo",nfile,fnm);
 489
 490         if ((ir->efep != efepNO || ir->bSimTemp) && ir->fepvals->nstdhdl > 0 &&
 491             (ir->fepvals->separate_dhdl_file == esepdhdlfileYES ) &&
 492             EI_DYNAMICS(ir->eI))
 493         {
 494             if (bAppendFiles)
 495             {
 496                 of->fp_dhdl = gmx_fio_fopen(opt2fn("-dhdl",nfile,fnm),filemode);
 497             }
 498             else
 499             {
 500                 of->fp_dhdl = open_dhdl(opt2fn("-dhdl",nfile,fnm),ir,oenv);
 501             }
 502         }
 503
 504         if (opt2bSet("-field",nfile,fnm) &&
 505             (ir->ex[XX].n || ir->ex[YY].n || ir->ex[ZZ].n))
 506         {
 507             if (bAppendFiles)
 508             {
 509                 of->fp_dhdl = gmx_fio_fopen(opt2fn("-field",nfile,fnm),
 510                                             filemode);
 511             }
 512             else
 513             {
 514                 of->fp_field = xvgropen(opt2fn("-field",nfile,fnm),
 515                                         "Applied electric field","Time (ps)",
 516                                         "E (V/nm)",oenv);
 517             }
 518         }
 519     }
 520
 521     return of;
 522 }
 523
 524 void done_mdoutf(gmx_mdoutf_t *of)
 525 {
 526     if (of->fp_ene != NULL)
 527     {
 528         close_enx(of->fp_ene);
 529     }
 530     if (of->fp_xtc)
 531     {
 532         close_xtc(of->fp_xtc);
 533     }
 534     if (of->fp_trn)
 535     {
 536         close_trn(of->fp_trn);
 537     }
 538     if (of->fp_dhdl != NULL)
 539     {
 540         gmx_fio_fclose(of->fp_dhdl);
 541     }
 542     if (of->fp_field != NULL)
 543     {
 544         gmx_fio_fclose(of->fp_field);
 545     }
 546
 547     sfree(of);
 548 }
 549
 550 void write_traj(FILE *fplog,t_commrec *cr,
 551                 gmx_mdoutf_t *of,
 552                 int mdof_flags,
 553                 gmx_mtop_t *top_global,
 554                 gmx_large_int_t step,double t,
 555                 t_state *state_local,t_state *state_global,
 556                 rvec *f_local,rvec *f_global,
 557                 int *n_xtc,rvec **x_xtc)
 558 {
 559     int     i,j;
 560     gmx_groups_t *groups;
 561     rvec    *xxtc;
 562     rvec *local_v;
 563     rvec *global_v;
 564
 565 #define MX(xvf) moveit(cr,GMX_LEFT,GMX_RIGHT,#xvf,xvf)
 566
 567     /* MRS -- defining these variables is to manage the difference
 568      * between half step and full step velocities, but there must be a better way . . . */
 569
 570     local_v  = state_local->v;
 571     global_v = state_global->v;
 572
 573     if (DOMAINDECOMP(cr))
 574     {
 575         if (mdof_flags & MDOF_CPT)
 576         {
 577             dd_collect_state(cr->dd,state_local,state_global);
 578         }
 579         else
 580         {
 581             if (mdof_flags & (MDOF_X | MDOF_XTC))
 582             {
 583                 dd_collect_vec(cr->dd,state_local,state_local->x,
 584                                state_global->x);
 585             }
 586             if (mdof_flags & MDOF_V)
 587             {
 588                 dd_collect_vec(cr->dd,state_local,local_v,
 589                                global_v);
 590             }
 591         }
 592         if (mdof_flags & MDOF_F)
 593         {
 594             dd_collect_vec(cr->dd,state_local,f_local,f_global);
 595         }
 596     }
 597     else
 598     {
 599         if (mdof_flags & MDOF_CPT)
 600         {
 601             /* All pointers in state_local are equal to state_global,
 602              * but we need to copy the non-pointer entries.
 603              */
 604             state_global->lambda = state_local->lambda;
 605             state_global->veta = state_local->veta;
 606             state_global->vol0 = state_local->vol0;
 607             copy_mat(state_local->box,state_global->box);
 608             copy_mat(state_local->boxv,state_global->boxv);
 609             copy_mat(state_local->svir_prev,state_global->svir_prev);
 610             copy_mat(state_local->fvir_prev,state_global->fvir_prev);
 611             copy_mat(state_local->pres_prev,state_global->pres_prev);
 612         }
 613         if (cr->nnodes > 1)
 614         {
 615             /* Particle decomposition, collect the data on the master node */
 616             if (mdof_flags & MDOF_CPT)
 617             {
 618                 if (state_local->flags & (1<<estX))   MX(state_global->x);
 619                 if (state_local->flags & (1<<estV))   MX(state_global->v);
 620                 if (state_local->flags & (1<<estSDX)) MX(state_global->sd_X);
 621                 if (state_global->nrngi > 1) {
 622                     if (state_local->flags & (1<<estLD_RNG)) {
 623 #ifdef GMX_MPI
 624                         MPI_Gather(state_local->ld_rng ,
 625                                    state_local->nrng*sizeof(state_local->ld_rng[0]),MPI_BYTE,
 626                                    state_global->ld_rng,
 627                                    state_local->nrng*sizeof(state_local->ld_rng[0]),MPI_BYTE,
 628                                    MASTERRANK(cr),cr->mpi_comm_mygroup);
 629 #endif
 630                     }
 631                     if (state_local->flags & (1<<estLD_RNGI))
 632                     {
 633 #ifdef GMX_MPI
 634                         MPI_Gather(state_local->ld_rngi,
 635                                    sizeof(state_local->ld_rngi[0]),MPI_BYTE,
 636                                    state_global->ld_rngi,
 637                                    sizeof(state_local->ld_rngi[0]),MPI_BYTE,
 638                                    MASTERRANK(cr),cr->mpi_comm_mygroup);
 639 #endif
 640                     }
 641                 }
 642             }
 643             else
 644             {
 645                 if (mdof_flags & (MDOF_X | MDOF_XTC)) MX(state_global->x);
 646                 if (mdof_flags & MDOF_V)              MX(global_v);
 647             }
 648             if (mdof_flags & MDOF_F) MX(f_global);
 649          }
 650      }
 651
 652      if (MASTER(cr))
 653      {
 654          if (mdof_flags & MDOF_CPT)
 655          {
 656              write_checkpoint(of->fn_cpt,of->bKeepAndNumCPT,
 657                               fplog,cr,of->eIntegrator,of->simulation_part,
 658                               of->bExpanded,of->elamstats,step,t,state_global);
 659          }
 660
 661          if (mdof_flags & (MDOF_X | MDOF_V | MDOF_F))
 662          {
 663             fwrite_trn(of->fp_trn,step,t,state_local->lambda[efptFEP],
 664                        state_local->box,top_global->natoms,
 665                        (mdof_flags & MDOF_X) ? state_global->x : NULL,
 666                        (mdof_flags & MDOF_V) ? global_v : NULL,
 667                        (mdof_flags & MDOF_F) ? f_global : NULL);
 668             if (gmx_fio_flush(of->fp_trn) != 0)
 669             {
 670                 gmx_file("Cannot write trajectory; maybe you are out of disk space?");
 671             }
 672             gmx_fio_check_file_position(of->fp_trn);
 673         }
 674         if (mdof_flags & MDOF_XTC) {
 675             groups = &top_global->groups;
 676             if (*n_xtc == -1)
 677             {
 678                 *n_xtc = 0;
 679                 for(i=0; (i<top_global->natoms); i++)
 680                 {
 681                     if (ggrpnr(groups,egcXTC,i) == 0)
 682                     {
 683                         (*n_xtc)++;
 684                     }
 685                 }
 686                 if (*n_xtc != top_global->natoms)
 687                 {
 688                     snew(*x_xtc,*n_xtc);
 689                 }
 690             }
 691             if (*n_xtc == top_global->natoms)
 692             {
 693                 xxtc = state_global->x;
 694             }
 695             else
 696             {
 697                 xxtc = *x_xtc;
 698                 j = 0;
 699                 for(i=0; (i<top_global->natoms); i++)
 700                 {
 701                     if (ggrpnr(groups,egcXTC,i) == 0)
 702                     {
 703                         copy_rvec(state_global->x[i],xxtc[j++]);
 704                     }
 705                 }
 706             }
 707             if (write_xtc(of->fp_xtc,*n_xtc,step,t,
 708                           state_local->box,xxtc,of->xtc_prec) == 0)
 709             {
 710                 gmx_fatal(FARGS,"XTC error - maybe you are out of disk space?");
 711             }
 712             gmx_fio_check_file_position(of->fp_xtc);
 713         }
 714     }
 715 }
 716