src/kernel/tpbconv.c

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