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