src/kernel/tpbconv.c

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