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