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 bool *bKeepIt(int gnx,int natoms,atom_id index[])
  62 {
  63   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(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[],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[],bool bKeep[],
 198                          t_ilist *il,int nratoms,const char *name)
 199 {
 200   t_iatom *ia;
 201   int i,j,newnr;
 202   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   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]1st.[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]2nd.[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     "Note 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]3rd.[bb] by creating a tpx file for a subset of your original",
 309     "tpx file, which is useful when you want to remove the solvent from",
 310     "your tpx file, or when you want to make e.g. a pure Ca tpx file.",
 311     "[BB]WARNING: this tpx file is not fully functional[bb].",
 312     "[BB]4th.[bb] by setting the charges of a specified group",
 313     "to zero. This is useful when doing free energy estimates",
 314     "using the LIE (Linear Interaction Energy) method."
 315   };
 316
 317   const char   *top_fn,*frame_fn;
 318   t_fileio     *fp;
 319   ener_file_t  fp_ener=NULL;
 320   t_trnheader head;
 321   int          i;
 322   gmx_large_int_t   nsteps_req,run_step,frame;
 323   double       run_t,state_t;
 324   bool         bOK,bNsteps,bExtend,bUntil,bTime,bTraj;
 325   bool         bFrame,bUse,bSel,bNeedEner,bReadEner,bScanEner;
 326   gmx_mtop_t   mtop;
 327   t_atoms      atoms;
 328   t_inputrec   *ir,*irnew=NULL;
 329   t_gromppopts *gopts;
 330   t_state      state;
 331   rvec         *newx=NULL,*newv=NULL,*tmpx,*tmpv;
 332   matrix       newbox;
 333   int          gnx;
 334   char         *grpname;
 335   atom_id      *index=NULL;
 336   int          nre;
 337   gmx_enxnm_t  *enm=NULL;
 338   t_enxframe   *fr_ener=NULL;
 339   char         buf[200],buf2[200];
 340   output_env_t oenv;
 341   t_filenm fnm[] = {
 342     { efTPX, NULL,  NULL,    ffREAD  },
 343     { efTRN, "-f",  NULL,    ffOPTRD },
 344     { efEDR, "-e",  NULL,    ffOPTRD },
 345     { efNDX, NULL,  NULL,    ffOPTRD },
 346     { efTPX, "-o",  "tpxout",ffWRITE }
 347   };
 348 #define NFILE asize(fnm)
 349
 350   /* Command line options */
 351   static int  nsteps_req_int = 0;
 352   static real start_t = -1.0, extend_t = 0.0, until_t = 0.0;
 353   static bool bContinuation = TRUE,bZeroQ = FALSE,bVel=TRUE;
 354   static t_pargs pa[] = {
 355     { "-extend",        FALSE, etREAL, {&extend_t},
 356       "Extend runtime by this amount (ps)" },
 357     { "-until",         FALSE, etREAL, {&until_t},
 358       "Extend runtime until this ending time (ps)" },
 359     { "-nsteps",        FALSE, etINT,  {&nsteps_req_int},
 360       "Change the number of steps" },
 361     { "-time",          FALSE, etREAL, {&start_t},
 362       "Continue from frame at this time (ps) instead of the last frame" },
 363     { "-zeroq",         FALSE, etBOOL, {&bZeroQ},
 364       "Set the charges of a group (from the index) to zero" },
 365     { "-vel",           FALSE, etBOOL, {&bVel},
 366       "Require velocities from trajectory" },
 367     { "-cont",          FALSE, etBOOL, {&bContinuation},
 368       "For exact continuation, the constraints should not be applied before the first step" }
 369   };
 370   int nerror = 0;
 371
 372   CopyRight(stdout,argv[0]);
 373
 374   /* Parse the command line */
 375   parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,
 376                     asize(desc),desc,0,NULL,&oenv);
 377
 378   /* Convert int to gmx_large_int_t */
 379   nsteps_req = nsteps_req_int;
 380   bNsteps = opt2parg_bSet("-nsteps",asize(pa),pa);
 381   bExtend = opt2parg_bSet("-extend",asize(pa),pa);
 382   bUntil  = opt2parg_bSet("-until",asize(pa),pa);
 383   bTime   = opt2parg_bSet("-time",asize(pa),pa);
 384   bTraj   = (opt2bSet("-f",NFILE,fnm) || bTime);
 385
 386   top_fn = ftp2fn(efTPX,NFILE,fnm);
 387   fprintf(stderr,"Reading toplogy and stuff from %s\n",top_fn);
 388
 389   snew(ir,1);
 390   read_tpx_state(top_fn,ir,&state,NULL,&mtop);
 391   run_step = ir->init_step;
 392   run_t    = ir->init_step*ir->delta_t + ir->init_t;
 393
 394   if (!EI_STATE_VELOCITY(ir->eI)) {
 395     bVel = FALSE;
 396   }
 397
 398   if (bTraj) {
 399     fprintf(stderr,"\n"
 400             "NOTE: Reading the state from trajectory is an obsolete feaure of tpbconv.\n"
 401             "      Continuation should be done by loading a checkpoint file with mdrun -cpi\n"
 402             "      This guarantees that all state variables are transferred.\n"
 403             "      tpbconv is now only useful for increasing nsteps,\n"
 404             "      but even that can often be avoided by using mdrun -maxh\n"
 405             "\n");
 406
 407     if (ir->bContinuation != bContinuation)
 408       fprintf(stderr,"Modifying ir->bContinuation to %s\n",
 409               bool_names[bContinuation]);
 410     ir->bContinuation = bContinuation;
 411
 412
 413     bNeedEner = (ir->epc == epcPARRINELLORAHMAN || ir->etc == etcNOSEHOOVER);
 414     bReadEner = (bNeedEner && ftp2bSet(efEDR,NFILE,fnm));
 415     bScanEner = (bReadEner && !bTime);
 416
 417     if (ir->epc != epcNO || EI_SD(ir->eI) || ir->eI == eiBD) {
 418       fprintf(stderr,"NOTE: The simulation uses pressure coupling and/or stochastic dynamics.\n"
 419               "tpbconv can not provide binary identical continuation.\n"
 420               "If you want that, supply a checkpoint file to mdrun\n\n");
 421     }
 422
 423     if (EI_SD(ir->eI) || ir->eI == eiBD) {
 424       fprintf(stderr,"\nChanging ld-seed from %d ",ir->ld_seed);
 425       ir->ld_seed = make_seed();
 426       fprintf(stderr,"to %d\n\n",ir->ld_seed);
 427     }
 428
 429     frame_fn = ftp2fn(efTRN,NFILE,fnm);
 430     fprintf(stderr,
 431             "\nREADING COORDS, VELS AND BOX FROM TRAJECTORY %s...\n\n",
 432             frame_fn);
 433
 434     fp = open_trn(frame_fn,"r");
 435     if (bScanEner) {
 436       fp_ener = open_enx(ftp2fn(efEDR,NFILE,fnm),"r");
 437       do_enxnms(fp_ener,&nre,&enm);
 438       snew(fr_ener,1);
 439       fr_ener->t = -1e-12;
 440     }
 441
 442     /* Now scan until the last set of x and v (step == 0)
 443      * or the ones at step step.
 444      */
 445     bFrame = TRUE;
 446     frame  = 0;
 447     while (bFrame) {
 448       bFrame = fread_trnheader(fp,&head,&bOK);
 449       if (bOK && frame == 0) {
 450         if (mtop.natoms != head.natoms)
 451           gmx_fatal(FARGS,"Number of atoms in Topology (%d) "
 452                       "is not the same as in Trajectory (%d)\n",
 453                       mtop.natoms,head.natoms);
 454         snew(newx,head.natoms);
 455         snew(newv,head.natoms);
 456       }
 457       bFrame = bFrame && bOK;
 458       if (bFrame) {
 459
 460         bOK = fread_htrn(fp,&head,newbox,newx,newv,NULL);
 461       }
 462       bFrame = bFrame && bOK;
 463       bUse = FALSE;
 464       if (bFrame &&
 465           (head.x_size) && (head.v_size || !bVel)) {
 466         bUse = TRUE;
 467         if (bScanEner) {
 468           /* Read until the energy time is >= the trajectory time */
 469           while (fr_ener->t < head.t && do_enx(fp_ener,fr_ener));
 470           bUse = (fr_ener->t == head.t);
 471         }
 472         if (bUse) {
 473           tmpx    = newx;
 474           newx    = state.x;
 475           state.x = tmpx;
 476           tmpv    = newv;
 477           newv    = state.v;
 478           state.v = tmpv;
 479           run_t        = head.t;
 480           run_step     = head.step;
 481           state.lambda = head.lambda;
 482           copy_mat(newbox,state.box);
 483         }
 484       }
 485       if (bFrame || !bOK) {
 486         sprintf(buf,"\r%s %s frame %s%s: step %s%s time %s",
 487                 "%s","%s","%6",gmx_large_int_fmt,"%6",gmx_large_int_fmt," %8.3f");
 488         fprintf(stderr,buf,
 489                 bUse ? "Read   " : "Skipped",ftp2ext(fn2ftp(frame_fn)),
 490                 frame,head.step,head.t);
 491         frame++;
 492         if (bTime && (head.t >= start_t))
 493           bFrame = FALSE;
 494       }
 495     }
 496     if (bScanEner) {
 497       close_enx(fp_ener);
 498       free_enxframe(fr_ener);
 499       free_enxnms(nre,enm);
 500     }
 501     close_trn(fp);
 502     fprintf(stderr,"\n");
 503
 504     if (!bOK)
 505       fprintf(stderr,"%s frame %s (step %s, time %g) is incomplete\n",
 506               ftp2ext(fn2ftp(frame_fn)),gmx_step_str(frame-1,buf2),
 507               gmx_step_str(head.step,buf),head.t);
 508     fprintf(stderr,"\nUsing frame of step %s time %g\n",
 509             gmx_step_str(run_step,buf),run_t);
 510
 511     if (bNeedEner) {
 512       if (bReadEner) {
 513         get_enx_state(ftp2fn(efEDR,NFILE,fnm),run_t,&mtop.groups,ir,&state);
 514       } else {
 515         fprintf(stderr,"\nWARNING: The simulation uses %s temperature and/or %s pressure coupling,\n"
 516                 "         the continuation will only be exact when an energy file is supplied\n\n",
 517                 ETCOUPLTYPE(etcNOSEHOOVER),
 518                 EPCOUPLTYPE(epcPARRINELLORAHMAN));
 519       }
 520     }
 521   }
 522
 523   if (bNsteps) {
 524     fprintf(stderr,"Setting nsteps to %s\n",gmx_step_str(nsteps_req,buf));
 525     ir->nsteps = nsteps_req;
 526   } else {
 527     /* Determine total number of steps remaining */
 528     if (bExtend) {
 529       ir->nsteps = ir->nsteps - (run_step - ir->init_step) + (int)(extend_t/ir->delta_t + 0.5);
 530       printf("Extending remaining runtime of by %g ps (now %s steps)\n",
 531              extend_t,gmx_step_str(ir->nsteps,buf));
 532     }
 533     else if (bUntil) {
 534       printf("nsteps = %s, run_step = %s, current_t = %g, until = %g\n",
 535              gmx_step_str(ir->nsteps,buf),
 536              gmx_step_str(run_step,buf2),
 537              run_t,until_t);
 538       ir->nsteps = (gmx_large_int_t)((until_t - run_t)/ir->delta_t + 0.5);
 539       printf("Extending remaining runtime until %g ps (now %s steps)\n",
 540              until_t,gmx_step_str(ir->nsteps,buf));
 541     }
 542     else {
 543       ir->nsteps -= run_step - ir->init_step;
 544       /* Print message */
 545       printf("%s steps (%g ps) remaining from first run.\n",
 546              gmx_step_str(ir->nsteps,buf),ir->nsteps*ir->delta_t);
 547
 548     }
 549   }
 550
 551   if (bNsteps || bZeroQ || (ir->nsteps > 0)) {
 552     ir->init_step = run_step;
 553
 554     if (ftp2bSet(efNDX,NFILE,fnm) ||
 555         !(bNsteps || bExtend || bUntil || bTraj)) {
 556       atoms = gmx_mtop_global_atoms(&mtop);
 557       get_index(&atoms,ftp2fn_null(efNDX,NFILE,fnm),1,
 558                 &gnx,&index,&grpname);
 559       if (!bZeroQ) {
 560         bSel = (gnx != state.natoms);
 561         for (i=0; ((i<gnx) && (!bSel)); i++)
 562           bSel = (i!=index[i]);
 563       }
 564       else
 565         bSel = FALSE;
 566       if (bSel) {
 567         fprintf(stderr,"Will write subset %s of original tpx containing %d "
 568                 "atoms\n",grpname,gnx);
 569         reduce_topology_x(gnx,index,&mtop,state.x,state.v);
 570         state.natoms = gnx;
 571       }
 572       else if (bZeroQ) {
 573         zeroq(gnx,index,&mtop);
 574         fprintf(stderr,"Zero-ing charges for group %s\n",grpname);
 575       }
 576       else
 577         fprintf(stderr,"Will write full tpx file (no selection)\n");
 578     }
 579
 580     state_t = ir->init_t + ir->init_step*ir->delta_t;
 581     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);
 582     fprintf(stderr,buf,ir->init_step,ir->nsteps);
 583     fprintf(stderr,"                                 time %10.3f and length %10.3f ps\n",
 584             state_t,ir->nsteps*ir->delta_t);
 585     write_tpx_state(opt2fn("-o",NFILE,fnm),ir,&state,&mtop);
 586   }
 587   else
 588     printf("You've simulated long enough. Not writing tpr file\n");
 589
 590   thanx(stderr);
 591
 592   return 0;
 593 }