src/programs/tpbconv/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 cmain (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,bFepState;
 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 int init_fep_state = 0;
 358   static gmx_bool bContinuation = TRUE,bZeroQ = FALSE,bVel=TRUE;
 359   static t_pargs pa[] = {
 360     { "-extend",        FALSE, etREAL, {&extend_t},
 361       "Extend runtime by this amount (ps)" },
 362     { "-until",         FALSE, etREAL, {&until_t},
 363       "Extend runtime until this ending time (ps)" },
 364     { "-nsteps",        FALSE, etINT,  {&nsteps_req_int},
 365       "Change the number of steps" },
 366     { "-time",          FALSE, etREAL, {&start_t},
 367       "Continue from frame at this time (ps) instead of the last frame" },
 368     { "-zeroq",         FALSE, etBOOL, {&bZeroQ},
 369       "Set the charges of a group (from the index) to zero" },
 370     { "-vel",           FALSE, etBOOL, {&bVel},
 371       "Require velocities from trajectory" },
 372     { "-cont",          FALSE, etBOOL, {&bContinuation},
 373       "For exact continuation, the constraints should not be applied before the first step" },
 374     { "-init_fep_state",FALSE, etINT, {&init_fep_state},
 375       "fep state to initialize from" },
 376   };
 377   int nerror = 0;
 378
 379   CopyRight(stderr,argv[0]);
 380
 381   /* Parse the command line */
 382   parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,
 383                     asize(desc),desc,0,NULL,&oenv);
 384
 385   /* Convert int to gmx_large_int_t */
 386   nsteps_req = nsteps_req_int;
 387   bNsteps   = opt2parg_bSet("-nsteps",asize(pa),pa);
 388   bExtend   = opt2parg_bSet("-extend",asize(pa),pa);
 389   bUntil    = opt2parg_bSet("-until",asize(pa),pa);
 390   bFepState = opt2parg_bSet("-init_fep_state",asize(pa),pa);
 391   bTime     = opt2parg_bSet("-time",asize(pa),pa);
 392   bTraj     = (opt2bSet("-f",NFILE,fnm) || bTime);
 393
 394   top_fn = ftp2fn(efTPX,NFILE,fnm);
 395   fprintf(stderr,"Reading toplogy and stuff from %s\n",top_fn);
 396
 397   snew(ir,1);
 398   read_tpx_state(top_fn,ir,&state,NULL,&mtop);
 399   run_step = ir->init_step;
 400   run_t    = ir->init_step*ir->delta_t + ir->init_t;
 401
 402   if (!EI_STATE_VELOCITY(ir->eI)) {
 403     bVel = FALSE;
 404   }
 405
 406   if (bTraj) {
 407     fprintf(stderr,"\n"
 408             "NOTE: Reading the state from trajectory is an obsolete feature of tpbconv.\n"
 409             "      Continuation should be done by loading a checkpoint file with mdrun -cpi\n"
 410             "      This guarantees that all state variables are transferred.\n"
 411             "      tpbconv is now only useful for increasing nsteps,\n"
 412             "      but even that can often be avoided by using mdrun -maxh\n"
 413             "\n");
 414
 415     if (ir->bContinuation != bContinuation)
 416       fprintf(stderr,"Modifying ir->bContinuation to %s\n",
 417               bool_names[bContinuation]);
 418     ir->bContinuation = bContinuation;
 419
 420
 421     bNeedEner = (ir->epc == epcPARRINELLORAHMAN || ir->etc == etcNOSEHOOVER);
 422     bReadEner = (bNeedEner && ftp2bSet(efEDR,NFILE,fnm));
 423     bScanEner = (bReadEner && !bTime);
 424
 425     if (ir->epc != epcNO || EI_SD(ir->eI) || ir->eI == eiBD) {
 426       fprintf(stderr,"NOTE: The simulation uses pressure coupling and/or stochastic dynamics.\n"
 427               "tpbconv can not provide binary identical continuation.\n"
 428               "If you want that, supply a checkpoint file to mdrun\n\n");
 429     }
 430
 431     if (EI_SD(ir->eI) || ir->eI == eiBD) {
 432       fprintf(stderr,"\nChanging ld-seed from %d ",ir->ld_seed);
 433       ir->ld_seed = make_seed();
 434       fprintf(stderr,"to %d\n\n",ir->ld_seed);
 435     }
 436
 437     frame_fn = ftp2fn(efTRN,NFILE,fnm);
 438
 439     if (fn2ftp(frame_fn) == efCPT)
 440     {
 441         int sim_part;
 442
 443         fprintf(stderr,
 444                 "\nREADING STATE FROM CHECKPOINT %s...\n\n",
 445                 frame_fn);
 446
 447         read_checkpoint_state(frame_fn,&sim_part,
 448                               &run_step,&run_t,&state);
 449     }
 450     else
 451     {
 452         fprintf(stderr,
 453                 "\nREADING COORDS, VELS AND BOX FROM TRAJECTORY %s...\n\n",
 454                 frame_fn);
 455
 456         fp = open_trn(frame_fn,"r");
 457         if (bScanEner)
 458         {
 459             fp_ener = open_enx(ftp2fn(efEDR,NFILE,fnm),"r");
 460             do_enxnms(fp_ener,&nre,&enm);
 461             snew(fr_ener,1);
 462             fr_ener->t = -1e-12;
 463         }
 464
 465         /* Now scan until the last set of x and v (step == 0)
 466          * or the ones at step step.
 467          */
 468         bFrame = TRUE;
 469         frame  = 0;
 470         while (bFrame)
 471         {
 472             bFrame = fread_trnheader(fp,&head,&bOK);
 473             if (bOK && frame == 0)
 474             {
 475                 if (mtop.natoms != head.natoms)
 476                     gmx_fatal(FARGS,"Number of atoms in Topology (%d) "
 477                               "is not the same as in Trajectory (%d)\n",
 478                               mtop.natoms,head.natoms);
 479                 snew(newx,head.natoms);
 480                 snew(newv,head.natoms);
 481             }
 482             bFrame = bFrame && bOK;
 483             if (bFrame)
 484             {
 485                 bOK = fread_htrn(fp,&head,newbox,newx,newv,NULL);
 486             }
 487             bFrame = bFrame && bOK;
 488             bUse = FALSE;
 489             if (bFrame &&
 490                 (head.x_size) && (head.v_size || !bVel))
 491             {
 492                 bUse = TRUE;
 493                 if (bScanEner)
 494                 {
 495                     /* Read until the energy time is >= the trajectory time */
 496                     while (fr_ener->t < head.t && do_enx(fp_ener,fr_ener));
 497                     bUse = (fr_ener->t == head.t);
 498                 }
 499                 if (bUse)
 500                 {
 501                     tmpx    = newx;
 502                     newx    = state.x;
 503                     state.x = tmpx;
 504                     tmpv    = newv;
 505                     newv    = state.v;
 506                     state.v = tmpv;
 507                     run_t        = head.t;
 508                     run_step     = head.step;
 509                     state.fep_state = head.fep_state;
 510                     state.lambda[efptFEP] = head.lambda;
 511                     copy_mat(newbox,state.box);
 512                 }
 513             }
 514             if (bFrame || !bOK)
 515             {
 516                 sprintf(buf,"\r%s %s frame %s%s: step %s%s time %s",
 517                         "%s","%s","%6",gmx_large_int_fmt,"%6",gmx_large_int_fmt," %8.3f");
 518                 fprintf(stderr,buf,
 519                         bUse ? "Read   " : "Skipped",ftp2ext(fn2ftp(frame_fn)),
 520                         frame,head.step,head.t);
 521                 frame++;
 522                 if (bTime && (head.t >= start_t))
 523                     bFrame = FALSE;
 524             }
 525         }
 526         if (bScanEner)
 527         {
 528             close_enx(fp_ener);
 529             free_enxframe(fr_ener);
 530             free_enxnms(nre,enm);
 531         }
 532         close_trn(fp);
 533         fprintf(stderr,"\n");
 534
 535         if (!bOK)
 536         {
 537             fprintf(stderr,"%s frame %s (step %s, time %g) is incomplete\n",
 538                     ftp2ext(fn2ftp(frame_fn)),gmx_step_str(frame-1,buf2),
 539                     gmx_step_str(head.step,buf),head.t);
 540         }
 541         fprintf(stderr,"\nUsing frame of step %s time %g\n",
 542                 gmx_step_str(run_step,buf),run_t);
 543
 544         if (bNeedEner)
 545         {
 546             if (bReadEner)
 547             {
 548                 get_enx_state(ftp2fn(efEDR,NFILE,fnm),run_t,&mtop.groups,ir,&state);
 549             }
 550             else
 551             {
 552                 fprintf(stderr,"\nWARNING: The simulation uses %s temperature and/or %s pressure coupling,\n"
 553                         "         the continuation will only be exact when an energy file is supplied\n\n",
 554                         ETCOUPLTYPE(etcNOSEHOOVER),
 555                         EPCOUPLTYPE(epcPARRINELLORAHMAN));
 556             }
 557         }
 558         if (bFepState)
 559         {
 560             ir->fepvals->init_fep_state = init_fep_state;
 561         }
 562     }
 563   }
 564
 565   if (bNsteps) {
 566     fprintf(stderr,"Setting nsteps to %s\n",gmx_step_str(nsteps_req,buf));
 567     ir->nsteps = nsteps_req;
 568   } else {
 569       /* Determine total number of steps remaining */
 570       if (bExtend)
 571       {
 572           ir->nsteps = ir->nsteps - (run_step - ir->init_step) + (gmx_large_int_t)(extend_t/ir->delta_t + 0.5);
 573           printf("Extending remaining runtime of by %g ps (now %s steps)\n",
 574                  extend_t,gmx_step_str(ir->nsteps,buf));
 575       }
 576       else if (bUntil)
 577       {
 578           printf("nsteps = %s, run_step = %s, current_t = %g, until = %g\n",
 579                  gmx_step_str(ir->nsteps,buf),
 580                  gmx_step_str(run_step,buf2),
 581                  run_t,until_t);
 582           ir->nsteps = (gmx_large_int_t)((until_t - run_t)/ir->delta_t + 0.5);
 583           printf("Extending remaining runtime until %g ps (now %s steps)\n",
 584                  until_t,gmx_step_str(ir->nsteps,buf));
 585       }
 586       else
 587       {
 588           ir->nsteps -= run_step - ir->init_step;
 589           /* Print message */
 590           printf("%s steps (%g ps) remaining from first run.\n",
 591                  gmx_step_str(ir->nsteps,buf),ir->nsteps*ir->delta_t);
 592       }
 593   }
 594
 595   if (bNsteps || bZeroQ || (ir->nsteps > 0))
 596   {
 597       ir->init_step = run_step;
 598
 599     if (ftp2bSet(efNDX,NFILE,fnm) ||
 600         !(bNsteps || bExtend || bUntil || bTraj))
 601     {
 602         atoms = gmx_mtop_global_atoms(&mtop);
 603         get_index(&atoms,ftp2fn_null(efNDX,NFILE,fnm),1,
 604                   &gnx,&index,&grpname);
 605         if (!bZeroQ)
 606         {
 607             bSel = (gnx != state.natoms);
 608             for (i=0; ((i<gnx) && (!bSel)); i++)
 609             {
 610                 bSel = (i!=index[i]);
 611             }
 612         }
 613         else
 614         {
 615             bSel = FALSE;
 616         }
 617         if (bSel)
 618         {
 619             fprintf(stderr,"Will write subset %s of original tpx containing %d "
 620                     "atoms\n",grpname,gnx);
 621             reduce_topology_x(gnx,index,&mtop,state.x,state.v);
 622             state.natoms = gnx;
 623         }
 624         else if (bZeroQ)
 625         {
 626             zeroq(gnx,index,&mtop);
 627             fprintf(stderr,"Zero-ing charges for group %s\n",grpname);
 628         }
 629         else
 630         {
 631             fprintf(stderr,"Will write full tpx file (no selection)\n");
 632         }
 633     }
 634
 635     state_t = ir->init_t + ir->init_step*ir->delta_t;
 636     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);
 637     fprintf(stderr,buf,ir->init_step,ir->nsteps);
 638     fprintf(stderr,"                                 time %10.3f and length %10.3f ps\n",
 639             state_t,ir->nsteps*ir->delta_t);
 640     write_tpx_state(opt2fn("-o",NFILE,fnm),ir,&state,&mtop);
 641   }
 642   else
 643   {
 644       printf("You've simulated long enough. Not writing tpr file\n");
 645   }
 646   thanx(stderr);
 647
 648   return 0;
 649 }