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             if (bPres)
 131             {
 132                 ato[to++] = afrom[from++];
 133             }
 134             break;
 135         default:
 136             if (bEner)
 137             {
 138                 ato[to++] = afrom[from++];
 139             }
 140             break;
 141         }
 142     }
 143
 144     return to;
 145 }
 146
 147 void global_stat(FILE *fplog,gmx_global_stat_t gs,
 148                  t_commrec *cr,gmx_enerdata_t *enerd,
 149                  tensor fvir,tensor svir,rvec mu_tot,
 150                  t_inputrec *inputrec,
 151                  gmx_ekindata_t *ekind,gmx_constr_t constr,
 152                  t_vcm *vcm,
 153                  int nsig,real *sig,
 154                  gmx_mtop_t *top_global, t_state *state_local,
 155                  gmx_bool bSumEkinhOld, int flags)
 156 /* instead of current system, gmx_booleans for summing virial, kinetic energy, and other terms */
 157 {
 158   t_bin  *rb;
 159   int    *itc0,*itc1;
 160   int    ie=0,ifv=0,isv=0,irmsd=0,imu=0;
 161   int    idedl=0,idvdll=0,idvdlnl=0,iepl=0,icm=0,imass=0,ica=0,inb=0;
 162   int    isig=-1;
 163   int    icj=-1,ici=-1,icx=-1;
 164   int    inn[egNR];
 165   real   copyenerd[F_NRE];
 166   int    nener,j;
 167   real   *rmsd_data=NULL;
 168   double nb;
 169   gmx_bool   bVV,bTemp,bEner,bPres,bConstrVir,bEkinAveVel,bFirstIterate,bReadEkin;
 170
 171   bVV           = EI_VV(inputrec->eI);
 172   bTemp         = flags & CGLO_TEMPERATURE;
 173   bEner         = flags & CGLO_ENERGY;
 174   bPres         = (flags & CGLO_PRESSURE);
 175   bConstrVir    = (flags & CGLO_CONSTRAINT);
 176   bFirstIterate = (flags & CGLO_FIRSTITERATE);
 177   bEkinAveVel   = (inputrec->eI==eiVV || (inputrec->eI==eiVVAK && bPres));
 178   bReadEkin     = (flags & CGLO_READEKIN);
 179
 180   rb   = gs->rb;
 181   itc0 = gs->itc0;
 182   itc1 = gs->itc1;
 183
 184
 185   reset_bin(rb);
 186   /* This routine copies all the data to be summed to one big buffer
 187    * using the t_bin struct.
 188    */
 189
 190   /* First, we neeed to identify which enerd->term should be
 191      communicated.  Temperature and pressure terms should only be
 192      communicated and summed when they need to be, to avoid repeating
 193      the sums and overcounting. */
 194
 195   nener = filter_enerdterm(enerd->term,TRUE,copyenerd,bTemp,bPres,bEner);
 196
 197   /* First, the data that needs to be communicated with velocity verlet every time
 198      This is just the constraint virial.*/
 199   if (bConstrVir) {
 200       isv = add_binr(rb,DIM*DIM,svir[0]);
 201       where();
 202   }
 203
 204 /* We need the force virial and the kinetic energy for the first time through with velocity verlet */
 205   if (bTemp || !bVV)
 206   {
 207       if (ekind)
 208       {
 209           for(j=0; (j<inputrec->opts.ngtc); j++)
 210           {
 211               if (bSumEkinhOld)
 212               {
 213                   itc0[j]=add_binr(rb,DIM*DIM,ekind->tcstat[j].ekinh_old[0]);
 214               }
 215               if (bEkinAveVel && !bReadEkin)
 216               {
 217                   itc1[j]=add_binr(rb,DIM*DIM,ekind->tcstat[j].ekinf[0]);
 218               }
 219               else if (!bReadEkin)
 220               {
 221                   itc1[j]=add_binr(rb,DIM*DIM,ekind->tcstat[j].ekinh[0]);
 222               }
 223           }
 224           /* these probably need to be put into one of these categories */
 225           where();
 226           idedl = add_binr(rb,1,&(ekind->dekindl));
 227           where();
 228           ica   = add_binr(rb,1,&(ekind->cosacc.mvcos));
 229           where();
 230       }
 231   }
 232   where();
 233
 234   if ((bPres || !bVV) && bFirstIterate)
 235   {
 236       ifv = add_binr(rb,DIM*DIM,fvir[0]);
 237   }
 238
 239
 240   if (bEner)
 241   {
 242       where();
 243       if (bFirstIterate)
 244       {
 245           ie  = add_binr(rb,nener,copyenerd);
 246       }
 247       where();
 248       if (constr)
 249       {
 250           rmsd_data = constr_rmsd_data(constr);
 251           if (rmsd_data)
 252           {
 253               irmsd = add_binr(rb,inputrec->eI==eiSD2 ? 3 : 2,rmsd_data);
 254           }
 255       }
 256       if (!NEED_MUTOT(*inputrec))
 257       {
 258           imu = add_binr(rb,DIM,mu_tot);
 259           where();
 260       }
 261
 262       if (bFirstIterate)
 263       {
 264           for(j=0; (j<egNR); j++)
 265           {
 266               inn[j]=add_binr(rb,enerd->grpp.nener,enerd->grpp.ener[j]);
 267           }
 268           where();
 269           if (inputrec->efep != efepNO)
 270           {
 271               idvdll  = add_bind(rb,efptNR,enerd->dvdl_lin);
 272               idvdlnl = add_bind(rb,efptNR,enerd->dvdl_nonlin);
 273               if (enerd->n_lambda > 0)
 274               {
 275                   iepl = add_bind(rb,enerd->n_lambda,enerd->enerpart_lambda);
 276               }
 277           }
 278       }
 279   }
 280
 281   if (vcm)
 282   {
 283       icm   = add_binr(rb,DIM*vcm->nr,vcm->group_p[0]);
 284       where();
 285       imass = add_binr(rb,vcm->nr,vcm->group_mass);
 286       where();
 287       if (vcm->mode == ecmANGULAR)
 288       {
 289           icj   = add_binr(rb,DIM*vcm->nr,vcm->group_j[0]);
 290           where();
 291           icx   = add_binr(rb,DIM*vcm->nr,vcm->group_x[0]);
 292           where();
 293           ici   = add_binr(rb,DIM*DIM*vcm->nr,vcm->group_i[0][0]);
 294           where();
 295       }
 296   }
 297
 298   if (DOMAINDECOMP(cr))
 299   {
 300       nb = cr->dd->nbonded_local;
 301       inb = add_bind(rb,1,&nb);
 302       }
 303   where();
 304   if (nsig > 0)
 305   {
 306       isig = add_binr(rb,nsig,sig);
 307   }
 308
 309   /* Global sum it all */
 310   if (debug)
 311   {
 312       fprintf(debug,"Summing %d energies\n",rb->maxreal);
 313   }
 314   sum_bin(rb,cr);
 315   where();
 316
 317   /* Extract all the data locally */
 318
 319   if (bConstrVir)
 320   {
 321       extract_binr(rb,isv ,DIM*DIM,svir[0]);
 322   }
 323
 324   /* We need the force virial and the kinetic energy for the first time through with velocity verlet */
 325   if (bTemp || !bVV)
 326   {
 327       if (ekind)
 328       {
 329           for(j=0; (j<inputrec->opts.ngtc); j++)
 330           {
 331               if (bSumEkinhOld)
 332               {
 333                   extract_binr(rb,itc0[j],DIM*DIM,ekind->tcstat[j].ekinh_old[0]);
 334               }
 335               if (bEkinAveVel && !bReadEkin) {
 336                   extract_binr(rb,itc1[j],DIM*DIM,ekind->tcstat[j].ekinf[0]);
 337               }
 338               else if (!bReadEkin)
 339               {
 340                   extract_binr(rb,itc1[j],DIM*DIM,ekind->tcstat[j].ekinh[0]);
 341               }
 342           }
 343           extract_binr(rb,idedl,1,&(ekind->dekindl));
 344           extract_binr(rb,ica,1,&(ekind->cosacc.mvcos));
 345           where();
 346       }
 347   }
 348   if ((bPres || !bVV) && bFirstIterate)
 349   {
 350       extract_binr(rb,ifv ,DIM*DIM,fvir[0]);
 351   }
 352
 353   if (bEner)
 354   {
 355       if (bFirstIterate)
 356       {
 357           extract_binr(rb,ie,nener,copyenerd);
 358           if (rmsd_data)
 359           {
 360               extract_binr(rb,irmsd,inputrec->eI==eiSD2 ? 3 : 2,rmsd_data);
 361           }
 362           if (!NEED_MUTOT(*inputrec))
 363           {
 364               extract_binr(rb,imu,DIM,mu_tot);
 365           }
 366
 367           for(j=0; (j<egNR); j++)
 368           {
 369               extract_binr(rb,inn[j],enerd->grpp.nener,enerd->grpp.ener[j]);
 370           }
 371           if (inputrec->efep != efepNO)
 372           {
 373               extract_bind(rb,idvdll ,efptNR,enerd->dvdl_lin);
 374               extract_bind(rb,idvdlnl,efptNR,enerd->dvdl_nonlin);
 375               if (enerd->n_lambda > 0)
 376               {
 377                   extract_bind(rb,iepl,enerd->n_lambda,enerd->enerpart_lambda);
 378               }
 379           }
 380           if (DOMAINDECOMP(cr))
 381           {
 382               extract_bind(rb,inb,1,&nb);
 383               if ((int)(nb + 0.5) != cr->dd->nbonded_global)
 384               {
 385                   dd_print_missing_interactions(fplog,cr,(int)(nb + 0.5),top_global,state_local);
 386               }
 387           }
 388           where();
 389
 390           filter_enerdterm(copyenerd,FALSE,enerd->term,bTemp,bPres,bEner);
 391       }
 392   }
 393
 394   if (vcm)
 395   {
 396       extract_binr(rb,icm,DIM*vcm->nr,vcm->group_p[0]);
 397       where();
 398       extract_binr(rb,imass,vcm->nr,vcm->group_mass);
 399       where();
 400       if (vcm->mode == ecmANGULAR)
 401       {
 402           extract_binr(rb,icj,DIM*vcm->nr,vcm->group_j[0]);
 403           where();
 404           extract_binr(rb,icx,DIM*vcm->nr,vcm->group_x[0]);
 405           where();
 406           extract_binr(rb,ici,DIM*DIM*vcm->nr,vcm->group_i[0][0]);
 407           where();
 408       }
 409   }
 410
 411   if (nsig > 0)
 412   {
 413       extract_binr(rb,isig,nsig,sig);
 414   }
 415   where();
 416 }
 417
 418 int do_per_step(gmx_large_int_t step,gmx_large_int_t nstep)
 419 {
 420   if (nstep != 0)
 421     return ((step % nstep)==0);
 422   else
 423     return 0;
 424 }
 425
 426 static void moveit(t_commrec *cr,
 427                    int left,int right,const char *s,rvec xx[])
 428 {
 429   if (!xx)
 430     return;
 431
 432   move_rvecs(cr,FALSE,FALSE,left,right,
 433              xx,NULL,(cr->nnodes-cr->npmenodes)-1,NULL);
 434 }
 435
 436 gmx_mdoutf_t *init_mdoutf(int nfile,const t_filenm fnm[],int mdrun_flags,
 437                           const t_commrec *cr,const t_inputrec *ir,
 438                           const output_env_t oenv)
 439 {
 440     gmx_mdoutf_t *of;
 441     char filemode[3];
 442     gmx_bool bAppendFiles;
 443
 444     snew(of,1);
 445
 446     of->fp_trn   = NULL;
 447     of->fp_ene   = NULL;
 448     of->fp_xtc   = NULL;
 449     of->fp_dhdl  = NULL;
 450     of->fp_field = NULL;
 451
 452     of->eIntegrator     = ir->eI;
 453     of->bExpanded       = ir->bExpanded;
 454     of->elamstats       = ir->expandedvals->elamstats;
 455     of->simulation_part = ir->simulation_part;
 456
 457     if (MASTER(cr))
 458     {
 459         bAppendFiles = (mdrun_flags & MD_APPENDFILES);
 460
 461         of->bKeepAndNumCPT = (mdrun_flags & MD_KEEPANDNUMCPT);
 462
 463         sprintf(filemode, bAppendFiles ? "a+" : "w+");
 464
 465         if ((EI_DYNAMICS(ir->eI) || EI_ENERGY_MINIMIZATION(ir->eI))
 466 #ifndef GMX_FAHCORE
 467             &&
 468             !(EI_DYNAMICS(ir->eI) &&
 469               ir->nstxout == 0 &&
 470               ir->nstvout == 0 &&
 471               ir->nstfout == 0)
 472 #endif
 473             )
 474         {
 475             of->fp_trn = open_trn(ftp2fn(efTRN,nfile,fnm), filemode);
 476         }
 477         if (EI_DYNAMICS(ir->eI) &&
 478             ir->nstxtcout > 0)
 479         {
 480             of->fp_xtc = open_xtc(ftp2fn(efXTC,nfile,fnm), filemode);
 481             of->xtc_prec = ir->xtcprec;
 482         }
 483         if (EI_DYNAMICS(ir->eI) || EI_ENERGY_MINIMIZATION(ir->eI))
 484         {
 485             of->fp_ene = open_enx(ftp2fn(efEDR,nfile,fnm), filemode);
 486         }
 487         of->fn_cpt = opt2fn("-cpo",nfile,fnm);
 488
 489         if ((ir->efep != efepNO || ir->bSimTemp) && ir->fepvals->nstdhdl > 0 &&
 490             (ir->fepvals->separate_dhdl_file == esepdhdlfileYES ) &&
 491             EI_DYNAMICS(ir->eI))
 492         {
 493             if (bAppendFiles)
 494             {
 495                 of->fp_dhdl = gmx_fio_fopen(opt2fn("-dhdl",nfile,fnm),filemode);
 496             }
 497             else
 498             {
 499                 of->fp_dhdl = open_dhdl(opt2fn("-dhdl",nfile,fnm),ir,oenv);
 500             }
 501         }
 502
 503         if (opt2bSet("-field",nfile,fnm) &&
 504             (ir->ex[XX].n || ir->ex[YY].n || ir->ex[ZZ].n))
 505         {
 506             if (bAppendFiles)
 507             {
 508                 of->fp_dhdl = gmx_fio_fopen(opt2fn("-field",nfile,fnm),
 509                                             filemode);
 510             }
 511             else
 512             {
 513                 of->fp_field = xvgropen(opt2fn("-field",nfile,fnm),
 514                                         "Applied electric field","Time (ps)",
 515                                         "E (V/nm)",oenv);
 516             }
 517         }
 518     }
 519
 520     return of;
 521 }
 522
 523 void done_mdoutf(gmx_mdoutf_t *of)
 524 {
 525     if (of->fp_ene != NULL)
 526     {
 527         close_enx(of->fp_ene);
 528     }
 529     if (of->fp_xtc)
 530     {
 531         close_xtc(of->fp_xtc);
 532     }
 533     if (of->fp_trn)
 534     {
 535         close_trn(of->fp_trn);
 536     }
 537     if (of->fp_dhdl != NULL)
 538     {
 539         gmx_fio_fclose(of->fp_dhdl);
 540     }
 541     if (of->fp_field != NULL)
 542     {
 543         gmx_fio_fclose(of->fp_field);
 544     }
 545
 546     sfree(of);
 547 }
 548
 549 void write_traj(FILE *fplog,t_commrec *cr,
 550                 gmx_mdoutf_t *of,
 551                 int mdof_flags,
 552                 gmx_mtop_t *top_global,
 553                 gmx_large_int_t step,double t,
 554                 t_state *state_local,t_state *state_global,
 555                 rvec *f_local,rvec *f_global,
 556                 int *n_xtc,rvec **x_xtc)
 557 {
 558     int     i,j;
 559     gmx_groups_t *groups;
 560     rvec    *xxtc;
 561     rvec *local_v;
 562     rvec *global_v;
 563
 564 #define MX(xvf) moveit(cr,GMX_LEFT,GMX_RIGHT,#xvf,xvf)
 565
 566     /* MRS -- defining these variables is to manage the difference
 567      * between half step and full step velocities, but there must be a better way . . . */
 568
 569     local_v  = state_local->v;
 570     global_v = state_global->v;
 571
 572     if (DOMAINDECOMP(cr))
 573     {
 574         if (mdof_flags & MDOF_CPT)
 575         {
 576             dd_collect_state(cr->dd,state_local,state_global);
 577         }
 578         else
 579         {
 580             if (mdof_flags & (MDOF_X | MDOF_XTC))
 581             {
 582                 dd_collect_vec(cr->dd,state_local,state_local->x,
 583                                state_global->x);
 584             }
 585             if (mdof_flags & MDOF_V)
 586             {
 587                 dd_collect_vec(cr->dd,state_local,local_v,
 588                                global_v);
 589             }
 590         }
 591         if (mdof_flags & MDOF_F)
 592         {
 593             dd_collect_vec(cr->dd,state_local,f_local,f_global);
 594         }
 595     }
 596     else
 597     {
 598         if (mdof_flags & MDOF_CPT)
 599         {
 600             /* All pointers in state_local are equal to state_global,
 601              * but we need to copy the non-pointer entries.
 602              */
 603             state_global->lambda = state_local->lambda;
 604             state_global->veta = state_local->veta;
 605             state_global->vol0 = state_local->vol0;
 606             copy_mat(state_local->box,state_global->box);
 607             copy_mat(state_local->boxv,state_global->boxv);
 608             copy_mat(state_local->svir_prev,state_global->svir_prev);
 609             copy_mat(state_local->fvir_prev,state_global->fvir_prev);
 610             copy_mat(state_local->pres_prev,state_global->pres_prev);
 611         }
 612         if (cr->nnodes > 1)
 613         {
 614             /* Particle decomposition, collect the data on the master node */
 615             if (mdof_flags & MDOF_CPT)
 616             {
 617                 if (state_local->flags & (1<<estX))   MX(state_global->x);
 618                 if (state_local->flags & (1<<estV))   MX(state_global->v);
 619                 if (state_local->flags & (1<<estSDX)) MX(state_global->sd_X);
 620                 if (state_global->nrngi > 1) {
 621                     if (state_local->flags & (1<<estLD_RNG)) {
 622 #ifdef GMX_MPI
 623                         MPI_Gather(state_local->ld_rng ,
 624                                    state_local->nrng*sizeof(state_local->ld_rng[0]),MPI_BYTE,
 625                                    state_global->ld_rng,
 626                                    state_local->nrng*sizeof(state_local->ld_rng[0]),MPI_BYTE,
 627                                    MASTERRANK(cr),cr->mpi_comm_mygroup);
 628 #endif
 629                     }
 630                     if (state_local->flags & (1<<estLD_RNGI))
 631                     {
 632 #ifdef GMX_MPI
 633                         MPI_Gather(state_local->ld_rngi,
 634                                    sizeof(state_local->ld_rngi[0]),MPI_BYTE,
 635                                    state_global->ld_rngi,
 636                                    sizeof(state_local->ld_rngi[0]),MPI_BYTE,
 637                                    MASTERRANK(cr),cr->mpi_comm_mygroup);
 638 #endif
 639                     }
 640                 }
 641             }
 642             else
 643             {
 644                 if (mdof_flags & (MDOF_X | MDOF_XTC)) MX(state_global->x);
 645                 if (mdof_flags & MDOF_V)              MX(global_v);
 646             }
 647             if (mdof_flags & MDOF_F) MX(f_global);
 648          }
 649      }
 650
 651      if (MASTER(cr))
 652      {
 653          if (mdof_flags & MDOF_CPT)
 654          {
 655              write_checkpoint(of->fn_cpt,of->bKeepAndNumCPT,
 656                               fplog,cr,of->eIntegrator,of->simulation_part,
 657                               of->bExpanded,of->elamstats,step,t,state_global);
 658          }
 659
 660          if (mdof_flags & (MDOF_X | MDOF_V | MDOF_F))
 661          {
 662             fwrite_trn(of->fp_trn,step,t,state_local->lambda[efptFEP],
 663                        state_local->box,top_global->natoms,
 664                        (mdof_flags & MDOF_X) ? state_global->x : NULL,
 665                        (mdof_flags & MDOF_V) ? global_v : NULL,
 666                        (mdof_flags & MDOF_F) ? f_global : NULL);
 667             if (gmx_fio_flush(of->fp_trn) != 0)
 668             {
 669                 gmx_file("Cannot write trajectory; maybe you are out of disk space?");
 670             }
 671             gmx_fio_check_file_position(of->fp_trn);
 672         }
 673         if (mdof_flags & MDOF_XTC) {
 674             groups = &top_global->groups;
 675             if (*n_xtc == -1)
 676             {
 677                 *n_xtc = 0;
 678                 for(i=0; (i<top_global->natoms); i++)
 679                 {
 680                     if (ggrpnr(groups,egcXTC,i) == 0)
 681                     {
 682                         (*n_xtc)++;
 683                     }
 684                 }
 685                 if (*n_xtc != top_global->natoms)
 686                 {
 687                     snew(*x_xtc,*n_xtc);
 688                 }
 689             }
 690             if (*n_xtc == top_global->natoms)
 691             {
 692                 xxtc = state_global->x;
 693             }
 694             else
 695             {
 696                 xxtc = *x_xtc;
 697                 j = 0;
 698                 for(i=0; (i<top_global->natoms); i++)
 699                 {
 700                     if (ggrpnr(groups,egcXTC,i) == 0)
 701                     {
 702                         copy_rvec(state_global->x[i],xxtc[j++]);
 703                     }
 704                 }
 705             }
 706             if (write_xtc(of->fp_xtc,*n_xtc,step,t,
 707                           state_local->box,xxtc,of->xtc_prec) == 0)
 708             {
 709                 gmx_fatal(FARGS,"XTC error - maybe you are out of disk space?");
 710             }
 711             gmx_fio_check_file_position(of->fp_xtc);
 712         }
 713     }
 714 }
 715