src/gromacs/mdlib/stat.c

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