src/kernel/tpbconv.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  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <math.h>
  41 #include "index.h"
  42 #include "gmx_fatal.h"
  43 #include "string2.h"
  44 #include "sysstuff.h"
  45 #include "smalloc.h"
  46 #include "macros.h"
  47 #include "names.h"
  48 #include "typedefs.h"
  49 #include "tpxio.h"
  50 #include "trnio.h"
  51 #include "enxio.h"
  52 #include "readir.h"
  53 #include "statutil.h"
  54 #include "copyrite.h"
  55 #include "futil.h"
  56 #include "vec.h"
  57 #include "mtop_util.h"
  58 #include "random.h"
  59 #include "checkpoint.h"
  60
  61 #define RANGECHK(i,n) if ((i)>=(n)) gmx_fatal(FARGS,"Your index file contains atomnumbers (e.g. %d)\nthat are larger than the number of atoms in the tpr file (%d)",(i),(n))
  62
  63 static gmx_bool *bKeepIt(int gnx,int natoms,atom_id index[])
  64 {
  65   gmx_bool *b;
  66   int  i;
  67
  68   snew(b,natoms);
  69   for(i=0; (i<gnx); i++) {
  70     RANGECHK(index[i],natoms);
  71     b[index[i]] = TRUE;
  72   }
  73
  74   return b;
  75 }
  76
  77 static atom_id *invind(int gnx,int natoms,atom_id index[])
  78 {
  79   atom_id *inv;
  80   int     i;
  81
  82   snew(inv,natoms);
  83   for(i=0; (i<gnx); i++) {
  84     RANGECHK(index[i],natoms);
  85     inv[index[i]] = i;
  86   }
  87
  88   return inv;
  89 }
  90
  91 static void reduce_block(gmx_bool bKeep[],t_block *block,
  92                          const char *name)
  93 {
  94   atom_id *index;
  95   int i,j,newi,newj;
  96
  97   snew(index,block->nr);
  98
  99   newi = newj = 0;
 100   for(i=0; (i<block->nr); i++) {
 101     for(j=block->index[i]; (j<block->index[i+1]); j++) {
 102       if (bKeep[j]) {
 103         newj++;
 104       }
 105     }
 106     if (newj > index[newi]) {
 107       newi++;
 108       index[newi] = newj;
 109     }
 110   }
 111
 112   fprintf(stderr,"Reduced block %8s from %6d to %6d index-, %6d to %6d a-entries\n",
 113           name,block->nr,newi,block->index[block->nr],newj);
 114   block->index = index;
 115   block->nr    = newi;
 116 }
 117
 118 static void reduce_blocka(atom_id invindex[],gmx_bool bKeep[],t_blocka *block,
 119                           const char *name)
 120 {
 121   atom_id *index,*a;
 122   int i,j,k,newi,newj;
 123
 124   snew(index,block->nr);
 125   snew(a,block->nra);
 126
 127   newi = newj = 0;
 128   for(i=0; (i<block->nr); i++) {
 129     for(j=block->index[i]; (j<block->index[i+1]); j++) {
 130       k = block->a[j];
 131       if (bKeep[k]) {
 132         a[newj] = invindex[k];
 133         newj++;
 134       }
 135     }
 136     if (newj > index[newi]) {
 137       newi++;
 138       index[newi] = newj;
 139     }
 140   }
 141
 142   fprintf(stderr,"Reduced block %8s from %6d to %6d index-, %6d to %6d a-entries\n",
 143           name,block->nr,newi,block->nra,newj);
 144   block->index = index;
 145   block->a     = a;
 146   block->nr    = newi;
 147   block->nra   = newj;
 148 }
 149
 150 static void reduce_rvec(int gnx,atom_id index[],rvec vv[])
 151 {
 152   rvec *ptr;
 153   int  i;
 154
 155   snew(ptr,gnx);
 156   for(i=0; (i<gnx); i++)
 157     copy_rvec(vv[index[i]],ptr[i]);
 158   for(i=0; (i<gnx); i++)
 159     copy_rvec(ptr[i],vv[i]);
 160   sfree(ptr);
 161 }
 162
 163 static void reduce_atom(int gnx,atom_id index[],t_atom atom[],char ***atomname,
 164                         int *nres,t_resinfo *resinfo)
 165 {
 166   t_atom *ptr;
 167   char   ***aname;
 168   t_resinfo *rinfo;
 169   int    i,nr;
 170
 171   snew(ptr,gnx);
 172   snew(aname,gnx);
 173   snew(rinfo,atom[index[gnx-1]].resind+1);
 174   for(i=0; (i<gnx); i++) {
 175     ptr[i]   = atom[index[i]];
 176     aname[i] = atomname[index[i]];
 177   }
 178   nr=-1;
 179   for(i=0; (i<gnx); i++) {
 180     atom[i]     = ptr[i];
 181     atomname[i] = aname[i];
 182     if ((i==0) || (atom[i].resind != atom[i-1].resind)) {
 183       nr++;
 184       rinfo[nr] = resinfo[atom[i].resind];
 185     }
 186     atom[i].resind = nr;
 187   }
 188   nr++;
 189   for(i=0; (i<nr); i++) {
 190     resinfo[i] = rinfo[i];
 191   }
 192   *nres=nr;
 193
 194   sfree(aname);
 195   sfree(ptr);
 196   sfree(rinfo);
 197 }
 198
 199 static void reduce_ilist(atom_id invindex[],gmx_bool bKeep[],
 200                          t_ilist *il,int nratoms,const char *name)
 201 {
 202   t_iatom *ia;
 203   int i,j,newnr;
 204   gmx_bool bB;
 205
 206   if (il->nr) {
 207     snew(ia,il->nr);
 208     newnr=0;
 209     for(i=0; (i<il->nr); i+=nratoms+1) {
 210       bB = TRUE;
 211       for(j=1; (j<=nratoms); j++) {
 212         bB = bB && bKeep[il->iatoms[i+j]];
 213       }
 214       if (bB) {
 215         ia[newnr++] = il->iatoms[i];
 216         for(j=1; (j<=nratoms); j++)
 217           ia[newnr++] = invindex[il->iatoms[i+j]];
 218       }
 219     }
 220     fprintf(stderr,"Reduced ilist %8s from %6d to %6d entries\n",
 221             name,il->nr/(nratoms+1),
 222           newnr/(nratoms+1));
 223
 224     il->nr = newnr;
 225     for(i=0; (i<newnr); i++)
 226       il->iatoms[i] = ia[i];
 227
 228     sfree(ia);
 229   }
 230 }
 231
 232 static void reduce_topology_x(int gnx,atom_id index[],
 233                               gmx_mtop_t *mtop,rvec x[],rvec v[])
 234 {
 235   t_topology top;
 236   gmx_bool    *bKeep;
 237   atom_id *invindex;
 238   int     i;
 239
 240   top = gmx_mtop_t_to_t_topology(mtop);
 241   bKeep    = bKeepIt(gnx,top.atoms.nr,index);
 242   invindex = invind(gnx,top.atoms.nr,index);
 243
 244   reduce_block(bKeep,&(top.cgs),"cgs");
 245   reduce_block(bKeep,&(top.mols),"mols");
 246   reduce_blocka(invindex,bKeep,&(top.excls),"excls");
 247   reduce_rvec(gnx,index,x);
 248   reduce_rvec(gnx,index,v);
 249   reduce_atom(gnx,index,top.atoms.atom,top.atoms.atomname,
 250               &(top.atoms.nres),top.atoms.resinfo);
 251
 252   for(i=0; (i<F_NRE); i++) {
 253     reduce_ilist(invindex,bKeep,&(top.idef.il[i]),
 254                  interaction_function[i].nratoms,
 255                  interaction_function[i].name);
 256   }
 257
 258   top.atoms.nr = gnx;
 259
 260   mtop->nmoltype = 1;
 261   snew(mtop->moltype,mtop->nmoltype);
 262   mtop->moltype[0].name  = mtop->name;
 263   mtop->moltype[0].atoms = top.atoms;
 264   for(i=0; i<F_NRE; i++) {
 265     mtop->moltype[0].ilist[i] = top.idef.il[i];
 266   }
 267   mtop->moltype[0].atoms = top.atoms;
 268   mtop->moltype[0].cgs   = top.cgs;
 269   mtop->moltype[0].excls = top.excls;
 270
 271   mtop->nmolblock = 1;
 272   snew(mtop->molblock,mtop->nmolblock);
 273   mtop->molblock[0].type       = 0;
 274   mtop->molblock[0].nmol       = 1;
 275   mtop->molblock[0].natoms_mol = top.atoms.nr;
 276   mtop->molblock[0].nposres_xA = 0;
 277   mtop->molblock[0].nposres_xB = 0;
 278
 279   mtop->natoms                 = top.atoms.nr;
 280 }
 281
 282 static void zeroq(int n,atom_id index[],gmx_mtop_t *mtop)
 283 {
 284   int mt,i;
 285
 286   for(mt=0; mt<mtop->nmoltype; mt++) {
 287     for(i=0; (i<mtop->moltype[mt].atoms.nr); i++) {
 288       mtop->moltype[mt].atoms.atom[index[i]].q = 0;
 289       mtop->moltype[mt].atoms.atom[index[i]].qB = 0;
 290     }
 291   }
 292 }
 293
 294 int main (int argc, char *argv[])
 295 {
 296   const char *desc[] = {
 297     "tpbconv can edit run input files in four ways.[PAR]",
 298     "[BB]1.[bb] by modifying the number of steps in a run input file",
 299     "with options [TT]-extend[tt], [TT]-until[tt] or [TT]-nsteps[tt]",
 300     "(nsteps=-1 means unlimited number of steps)[PAR]",
 301     "[BB]2.[bb] (OBSOLETE) by creating a run input file",
 302     "for a continuation run when your simulation has crashed due to e.g.",
 303     "a full disk, or by making a continuation run input file.",
 304     "This option is obsolete, since mdrun now writes and reads",
 305     "checkpoint files.",
 306     "[BB]Note[bb] that a frame with coordinates and velocities is needed.",
 307     "When pressure and/or Nose-Hoover temperature coupling is used",
 308     "an energy file can be supplied to get an exact continuation",
 309     "of the original run.[PAR]",
 310     "[BB]3.[bb] by creating a [TT].tpx[tt] file for a subset of your original",
 311     "tpx file, which is useful when you want to remove the solvent from",
 312     "your [TT].tpx[tt] file, or when you want to make e.g. a pure C[GRK]alpha[grk] [TT].tpx[tt] file.",
 313     "Note that you may need to use [TT]-nsteps -1[tt] (or similar) to get",
 314     "this to work.",
 315     "[BB]WARNING: this [TT].tpx[tt] file is not fully functional[bb].[PAR]",
 316     "[BB]4.[bb] by setting the charges of a specified group",
 317     "to zero. This is useful when doing free energy estimates",
 318     "using the LIE (Linear Interaction Energy) method."
 319   };
 320
 321   const char   *top_fn,*frame_fn;
 322   t_fileio     *fp;
 323   ener_file_t  fp_ener=NULL;
 324   t_trnheader head;
 325   int          i;
 326   gmx_large_int_t   nsteps_req,run_step,frame;
 327   double       run_t,state_t;
 328   gmx_bool         bOK,bNsteps,bExtend,bUntil,bTime,bTraj;
 329   gmx_bool         bFrame,bUse,bSel,bNeedEner,bReadEner,bScanEner;
 330   gmx_mtop_t   mtop;
 331   t_atoms      atoms;
 332   t_inputrec   *ir,*irnew=NULL;
 333   t_gromppopts *gopts;
 334   t_state      state;
 335   rvec         *newx=NULL,*newv=NULL,*tmpx,*tmpv;
 336   matrix       newbox;
 337   int          gnx;
 338   char         *grpname;
 339   atom_id      *index=NULL;
 340   int          nre;
 341   gmx_enxnm_t  *enm=NULL;
 342   t_enxframe   *fr_ener=NULL;
 343   char         buf[200],buf2[200];
 344   output_env_t oenv;
 345   t_filenm fnm[] = {
 346     { efTPX, NULL,  NULL,    ffREAD  },
 347     { efTRN, "-f",  NULL,    ffOPTRD },
 348     { efEDR, "-e",  NULL,    ffOPTRD },
 349     { efNDX, NULL,  NULL,    ffOPTRD },
 350     { efTPX, "-o",  "tpxout",ffWRITE }
 351   };
 352 #define NFILE asize(fnm)
 353
 354   /* Command line options */
 355   static int  nsteps_req_int = 0;
 356   static real start_t = -1.0, extend_t = 0.0, until_t = 0.0;
 357   static gmx_bool bContinuation = TRUE,bZeroQ = FALSE,bVel=TRUE;
 358   static t_pargs pa[] = {
 359     { "-extend",        FALSE, etREAL, {&extend_t},
 360       "Extend runtime by this amount (ps)" },
 361     { "-until",         FALSE, etREAL, {&until_t},
 362       "Extend runtime until this ending time (ps)" },
 363     { "-nsteps",        FALSE, etINT,  {&nsteps_req_int},
 364       "Change the number of steps" },
 365     { "-time",          FALSE, etREAL, {&start_t},
 366       "Continue from frame at this time (ps) instead of the last frame" },
 367     { "-zeroq",         FALSE, etBOOL, {&bZeroQ},
 368       "Set the charges of a group (from the index) to zero" },
 369     { "-vel",           FALSE, etBOOL, {&bVel},
 370       "Require velocities from trajectory" },
 371     { "-cont",          FALSE, etBOOL, {&bContinuation},
 372       "For exact continuation, the constraints should not be applied before the first step" }
 373   };
 374   int nerror = 0;
 375
 376   CopyRight(stdout,argv[0]);
 377
 378   /* Parse the command line */
 379   parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,
 380                     asize(desc),desc,0,NULL,&oenv);
 381
 382   /* Convert int to gmx_large_int_t */
 383   nsteps_req = nsteps_req_int;
 384   bNsteps = opt2parg_bSet("-nsteps",asize(pa),pa);
 385   bExtend = opt2parg_bSet("-extend",asize(pa),pa);
 386   bUntil  = opt2parg_bSet("-until",asize(pa),pa);
 387   bTime   = opt2parg_bSet("-time",asize(pa),pa);
 388   bTraj   = (opt2bSet("-f",NFILE,fnm) || bTime);
 389
 390   top_fn = ftp2fn(efTPX,NFILE,fnm);
 391   fprintf(stderr,"Reading toplogy and stuff from %s\n",top_fn);
 392
 393   snew(ir,1);
 394   read_tpx_state(top_fn,ir,&state,NULL,&mtop);
 395   run_step = ir->init_step;
 396   run_t    = ir->init_step*ir->delta_t + ir->init_t;
 397
 398   if (!EI_STATE_VELOCITY(ir->eI)) {
 399     bVel = FALSE;
 400   }
 401
 402   if (bTraj) {
 403     fprintf(stderr,"\n"
 404             "NOTE: Reading the state from trajectory is an obsolete feature of tpbconv.\n"
 405             "      Continuation should be done by loading a checkpoint file with mdrun -cpi\n"
 406             "      This guarantees that all state variables are transferred.\n"
 407             "      tpbconv is now only useful for increasing nsteps,\n"
 408             "      but even that can often be avoided by using mdrun -maxh\n"
 409             "\n");
 410
 411     if (ir->bContinuation != bContinuation)
 412       fprintf(stderr,"Modifying ir->bContinuation to %s\n",
 413               bool_names[bContinuation]);
 414     ir->bContinuation = bContinuation;
 415
 416
 417     bNeedEner = (ir->epc == epcPARRINELLORAHMAN || ir->etc == etcNOSEHOOVER);
 418     bReadEner = (bNeedEner && ftp2bSet(efEDR,NFILE,fnm));
 419     bScanEner = (bReadEner && !bTime);
 420
 421     if (ir->epc != epcNO || EI_SD(ir->eI) || ir->eI == eiBD) {
 422       fprintf(stderr,"NOTE: The simulation uses pressure coupling and/or stochastic dynamics.\n"
 423               "tpbconv can not provide binary identical continuation.\n"
 424               "If you want that, supply a checkpoint file to mdrun\n\n");
 425     }
 426
 427     if (EI_SD(ir->eI) || ir->eI == eiBD) {
 428       fprintf(stderr,"\nChanging ld-seed from %d ",ir->ld_seed);
 429       ir->ld_seed = make_seed();
 430       fprintf(stderr,"to %d\n\n",ir->ld_seed);
 431     }
 432
 433     frame_fn = ftp2fn(efTRN,NFILE,fnm);
 434
 435         if (fn2ftp(frame_fn) == efCPT)
 436         {
 437             int sim_part;
 438
 439             fprintf(stderr,
 440                     "\nREADING STATE FROM CHECKPOINT %s...\n\n",
 441                     frame_fn);
 442
 443             read_checkpoint_state(frame_fn,&sim_part,
 444                                   &run_step,&run_t,&state);
 445         }
 446         else
 447         {
 448             fprintf(stderr,
 449                     "\nREADING COORDS, VELS AND BOX FROM TRAJECTORY %s...\n\n",
 450                     frame_fn);
 451
 452             fp = open_trn(frame_fn,"r");
 453             if (bScanEner)
 454             {
 455                 fp_ener = open_enx(ftp2fn(efEDR,NFILE,fnm),"r");
 456                 do_enxnms(fp_ener,&nre,&enm);
 457                 snew(fr_ener,1);
 458                 fr_ener->t = -1e-12;
 459             }
 460
 461             /* Now scan until the last set of x and v (step == 0)
 462              * or the ones at step step.
 463              */
 464             bFrame = TRUE;
 465             frame  = 0;
 466             while (bFrame)
 467             {
 468                 bFrame = fread_trnheader(fp,&head,&bOK);
 469                 if (bOK && frame == 0)
 470                 {
 471                     if (mtop.natoms != head.natoms)
 472                         gmx_fatal(FARGS,"Number of atoms in Topology (%d) "
 473                                   "is not the same as in Trajectory (%d)\n",
 474                                   mtop.natoms,head.natoms);
 475                     snew(newx,head.natoms);
 476                     snew(newv,head.natoms);
 477                 }
 478                 bFrame = bFrame && bOK;
 479                 if (bFrame)
 480                 {
 481                     bOK = fread_htrn(fp,&head,newbox,newx,newv,NULL);
 482                 }
 483                 bFrame = bFrame && bOK;
 484                 bUse = FALSE;
 485                 if (bFrame &&
 486                     (head.x_size) && (head.v_size || !bVel))
 487                 {
 488                     bUse = TRUE;
 489                     if (bScanEner)
 490                     {
 491                         /* Read until the energy time is >= the trajectory time */
 492                         while (fr_ener->t < head.t && do_enx(fp_ener,fr_ener));
 493                         bUse = (fr_ener->t == head.t);
 494                     }
 495                     if (bUse)
 496                     {
 497                         tmpx    = newx;
 498                         newx    = state.x;
 499                         state.x = tmpx;
 500                         tmpv    = newv;
 501                         newv    = state.v;
 502                         state.v = tmpv;
 503                         run_t        = head.t;
 504                         run_step     = head.step;
 505                         state.lambda = head.lambda;
 506                         copy_mat(newbox,state.box);
 507                     }
 508                 }
 509                 if (bFrame || !bOK)
 510                 {
 511                     sprintf(buf,"\r%s %s frame %s%s: step %s%s time %s",
 512                             "%s","%s","%6",gmx_large_int_fmt,"%6",gmx_large_int_fmt," %8.3f");
 513                     fprintf(stderr,buf,
 514                             bUse ? "Read   " : "Skipped",ftp2ext(fn2ftp(frame_fn)),
 515                             frame,head.step,head.t);
 516                     frame++;
 517                     if (bTime && (head.t >= start_t))
 518                         bFrame = FALSE;
 519                 }
 520             }
 521             if (bScanEner)
 522             {
 523                 close_enx(fp_ener);
 524                 free_enxframe(fr_ener);
 525                 free_enxnms(nre,enm);
 526             }
 527             close_trn(fp);
 528             fprintf(stderr,"\n");
 529
 530             if (!bOK)
 531             {
 532                 fprintf(stderr,"%s frame %s (step %s, time %g) is incomplete\n",
 533                         ftp2ext(fn2ftp(frame_fn)),gmx_step_str(frame-1,buf2),
 534                         gmx_step_str(head.step,buf),head.t);
 535             }
 536             fprintf(stderr,"\nUsing frame of step %s time %g\n",
 537                     gmx_step_str(run_step,buf),run_t);
 538
 539             if (bNeedEner)
 540             {
 541                 if (bReadEner)
 542                 {
 543                     get_enx_state(ftp2fn(efEDR,NFILE,fnm),run_t,&mtop.groups,ir,&state);
 544                 }
 545                 else
 546                 {
 547                     fprintf(stderr,"\nWARNING: The simulation uses %s temperature and/or %s pressure coupling,\n"
 548                             "         the continuation will only be exact when an energy file is supplied\n\n",
 549                             ETCOUPLTYPE(etcNOSEHOOVER),
 550                             EPCOUPLTYPE(epcPARRINELLORAHMAN));
 551                 }
 552             }
 553         }
 554   }
 555
 556   if (bNsteps) {
 557     fprintf(stderr,"Setting nsteps to %s\n",gmx_step_str(nsteps_req,buf));
 558     ir->nsteps = nsteps_req;
 559   } else {
 560     /* Determine total number of steps remaining */
 561     if (bExtend) {
 562       ir->nsteps = ir->nsteps - (run_step - ir->init_step) + (gmx_large_int_t)(extend_t/ir->delta_t + 0.5);
 563       printf("Extending remaining runtime of by %g ps (now %s steps)\n",
 564              extend_t,gmx_step_str(ir->nsteps,buf));
 565     }
 566     else if (bUntil) {
 567       printf("nsteps = %s, run_step = %s, current_t = %g, until = %g\n",
 568              gmx_step_str(ir->nsteps,buf),
 569              gmx_step_str(run_step,buf2),
 570              run_t,until_t);
 571       ir->nsteps = (gmx_large_int_t)((until_t - run_t)/ir->delta_t + 0.5);
 572       printf("Extending remaining runtime until %g ps (now %s steps)\n",
 573              until_t,gmx_step_str(ir->nsteps,buf));
 574     }
 575     else {
 576       ir->nsteps -= run_step - ir->init_step;
 577       /* Print message */
 578       printf("%s steps (%g ps) remaining from first run.\n",
 579              gmx_step_str(ir->nsteps,buf),ir->nsteps*ir->delta_t);
 580
 581     }
 582   }
 583
 584   if (bNsteps || bZeroQ || (ir->nsteps > 0)) {
 585     ir->init_step = run_step;
 586
 587     if (ftp2bSet(efNDX,NFILE,fnm) ||
 588         !(bNsteps || bExtend || bUntil || bTraj)) {
 589       atoms = gmx_mtop_global_atoms(&mtop);
 590       get_index(&atoms,ftp2fn_null(efNDX,NFILE,fnm),1,
 591                 &gnx,&index,&grpname);
 592       if (!bZeroQ) {
 593         bSel = (gnx != state.natoms);
 594         for (i=0; ((i<gnx) && (!bSel)); i++)
 595           bSel = (i!=index[i]);
 596       }
 597       else
 598         bSel = FALSE;
 599       if (bSel) {
 600         fprintf(stderr,"Will write subset %s of original tpx containing %d "
 601                 "atoms\n",grpname,gnx);
 602         reduce_topology_x(gnx,index,&mtop,state.x,state.v);
 603         state.natoms = gnx;
 604       }
 605       else if (bZeroQ) {
 606         zeroq(gnx,index,&mtop);
 607         fprintf(stderr,"Zero-ing charges for group %s\n",grpname);
 608       }
 609       else
 610         fprintf(stderr,"Will write full tpx file (no selection)\n");
 611     }
 612
 613     state_t = ir->init_t + ir->init_step*ir->delta_t;
 614     sprintf(buf,   "Writing statusfile with starting step %s%s and length %s%s steps...\n","%10",gmx_large_int_fmt,"%10",gmx_large_int_fmt);
 615     fprintf(stderr,buf,ir->init_step,ir->nsteps);
 616     fprintf(stderr,"                                 time %10.3f and length %10.3f ps\n",
 617             state_t,ir->nsteps*ir->delta_t);
 618     write_tpx_state(opt2fn("-o",NFILE,fnm),ir,&state,&mtop);
 619   }
 620   else
 621     printf("You've simulated long enough. Not writing tpr file\n");
 622
 623   thanx(stderr);
 624
 625   return 0;
 626 }