src/kernel/tpbconv.c

   1 /*
   2  * $Id$
   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 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_tpbconv_c = "$Id$";
  30
  31 #include <math.h>
  32 #include "rdgroup.h"
  33 #include "fatal.h"
  34 #include "string2.h"
  35 #include "sysstuff.h"
  36 #include "smalloc.h"
  37 #include "macros.h"
  38 #include "names.h"
  39 #include "typedefs.h"
  40 #include "tpxio.h"
  41 #include "trnio.h"
  42 #include "readir.h"
  43 #include "statutil.h"
  44 #include "copyrite.h"
  45 #include "futil.h"
  46 #include "assert.h"
  47 #include "vec.h"
  48 #include "rdgroup.h"
  49
  50 static bool *bKeepIt(int gnx,int natoms,atom_id index[])
  51 {
  52   bool *b;
  53   int  i;
  54
  55   snew(b,natoms);
  56   for(i=0; (i<gnx); i++) {
  57     assert(index[i] < natoms);
  58     b[index[i]] = TRUE;
  59   }
  60
  61   return b;
  62 }
  63
  64 static atom_id *invind(int gnx,int natoms,atom_id index[])
  65 {
  66   atom_id *inv;
  67   int     i;
  68
  69   snew(inv,natoms);
  70   for(i=0; (i<gnx); i++) {
  71     assert(index[i] < natoms);
  72     inv[index[i]] = i;
  73   }
  74
  75   return inv;
  76 }
  77
  78 static void reduce_block(atom_id invindex[],bool bKeep[],t_block *block,
  79                          char *name,bool bExcl)
  80 {
  81   atom_id *index,*a;
  82   int i,j,k,newi,newj;
  83
  84   snew(index,block->nr);
  85   snew(a,block->nra);
  86
  87   newi = newj = 0;
  88   for(i=0; (i<block->nr); i++) {
  89     if (!bExcl || bKeep[i]) {
  90       for(j=block->index[i]; (j<block->index[i+1]); j++) {
  91         k = block->a[j];
  92         if (bKeep[k]) {
  93           a[newj] = invindex[k];
  94           newj++;
  95         }
  96       }
  97       if (newj > index[newi]) {
  98         newi++;
  99         index[newi] = newj;
 100       }
 101     }
 102   }
 103
 104   fprintf(stderr,"Reduced block %8s from %6d to %6d index-, %6d to %6d a-entries\n",
 105           name,block->nr,newi,block->nra,newj);
 106   block->index = index;
 107   block->a     = a;
 108   block->nr    = newi;
 109   block->nra   = newj;
 110 }
 111
 112 static void reduce_rvec(int gnx,atom_id index[],rvec vv[])
 113 {
 114   rvec *ptr;
 115   int  i;
 116
 117   snew(ptr,gnx);
 118   for(i=0; (i<gnx); i++)
 119     copy_rvec(vv[index[i]],ptr[i]);
 120   for(i=0; (i<gnx); i++)
 121     copy_rvec(ptr[i],vv[i]);
 122   sfree(ptr);
 123 }
 124
 125 static void reduce_atom(int gnx,atom_id index[],t_atom atom[],char ***atomname,
 126                         int *nres, char ***resname)
 127 {
 128   t_atom *ptr;
 129   char   ***aname,***rname;
 130   int    i,nr;
 131
 132   snew(ptr,gnx);
 133   snew(aname,gnx);
 134   snew(rname,atom[index[gnx-1]].resnr+1);
 135   for(i=0; (i<gnx); i++) {
 136     ptr[i]   = atom[index[i]];
 137     aname[i] = atomname[index[i]];
 138   }
 139   nr=-1;
 140   for(i=0; (i<gnx); i++) {
 141     atom[i]     = ptr[i];
 142     atomname[i] = aname[i];
 143     if ((i==0) || (atom[i].resnr != atom[i-1].resnr)) {
 144       nr++;
 145       rname[nr]=resname[atom[i].resnr];
 146     }
 147     atom[i].resnr=nr;
 148   }
 149   nr++;
 150   for(i=0; (i<nr); i++)
 151     resname[i]=rname[i];
 152   *nres=nr;
 153
 154   sfree(aname);
 155   sfree(ptr);
 156   sfree(rname);
 157 }
 158
 159 static void reduce_ilist(atom_id invindex[],bool bKeep[],
 160                          t_ilist *il,int nratoms,char *name)
 161 {
 162   t_iatom *ia;
 163   int i,j,newnr;
 164   bool bB;
 165
 166   if (il->nr) {
 167     snew(ia,il->nr);
 168     newnr=0;
 169     for(i=0; (i<il->nr); i+=nratoms+1) {
 170       bB = TRUE;
 171       for(j=1; (j<=nratoms); j++) {
 172         bB = bB && bKeep[il->iatoms[i+j]];
 173       }
 174       if (bB) {
 175         ia[newnr++] = il->iatoms[i];
 176         for(j=1; (j<=nratoms); j++)
 177           ia[newnr++] = invindex[il->iatoms[i+j]];
 178       }
 179     }
 180     fprintf(stderr,"Reduced ilist %8s from %6d to %6d entries\n",
 181             name,il->nr/(nratoms+1),
 182           newnr/(nratoms+1));
 183
 184     il->nr = newnr;
 185     for(i=0; (i<newnr); i++)
 186       il->iatoms[i] = ia[i];
 187
 188     for(i=0; (i<MAXPROC); i++)
 189       il->multinr[i] = newnr;
 190
 191     sfree(ia);
 192   }
 193 }
 194
 195 static void reduce_topology_x(int gnx,atom_id index[],
 196                               t_topology *top,rvec x[],rvec v[])
 197 {
 198   bool    *bKeep;
 199   atom_id *invindex;
 200   int     i;
 201
 202   bKeep    = bKeepIt(gnx,top->atoms.nr,index);
 203   invindex = invind(gnx,top->atoms.nr,index);
 204
 205   for(i=0; (i<ebNR); i++)
 206     reduce_block(invindex,bKeep,&(top->blocks[i]),eblock_names[i],FALSE);
 207   reduce_block(invindex,bKeep,&(top->atoms.excl),"EXCL",TRUE);
 208   reduce_rvec(gnx,index,x);
 209   reduce_rvec(gnx,index,v);
 210   reduce_atom(gnx,index,top->atoms.atom,top->atoms.atomname,
 211               &(top->atoms.nres),top->atoms.resname);
 212
 213   for(i=0; (i<F_NRE); i++) {
 214     reduce_ilist(invindex,bKeep,&(top->idef.il[i]),
 215                  interaction_function[i].nratoms,
 216                  interaction_function[i].name);
 217   }
 218
 219   top->atoms.nr = gnx;
 220 }
 221
 222 int main (int argc, char *argv[])
 223 {
 224   static char *desc[] = {
 225     "tpbconv can edit run input files in two ways.[PAR]"
 226     "[BB]1st.[bb] by creating a run input file",
 227     "for a continuation run when your simulation has crashed due to e.g.",
 228     "a full disk, or by making a continuation run input file.",
 229     "Note that a frame with coordinates and velocities is needed,",
 230     "which means that when you never write velocities, you can not use",
 231     "tpbconv and you have to start the run again from the beginning.[PAR]",
 232     "[BB]2nd.[bb] by creating a tpx file for a subset of your original",
 233     "tpx file, which is useful when you want to remove the solvent from",
 234     "your tpx file, or when you want to make e.g. a pure Ca tpx file.",
 235     "[BB]WARNING: this tpx file is not fully functional[bb]."
 236   };
 237
 238   char         *top_fn,*frame_fn;
 239   int          fp;
 240   t_tpxheader  tpx;
 241   t_trnheader head;
 242   int          i,natoms,frame,step,run_step;
 243   real         run_t,run_lambda;
 244   bool         bOK,bFrame,bTime;
 245   t_topology   top;
 246   t_inputrec   *ir,*irnew=NULL;
 247   t_gromppopts *gopts;
 248   rvec         *x=NULL,*v=NULL,*newx,*newv,*tmpx,*tmpv;
 249   matrix       box,newbox;
 250   int          gnx;
 251   char         *grpname;
 252   atom_id      *index=NULL;
 253   t_filenm fnm[] = {
 254     { efTPX, NULL,  NULL,    ffREAD  },
 255     { efTRN, "-f",  NULL,    ffOPTRD },
 256     { efNDX, NULL,  NULL,    ffOPTRD },
 257     { efTPX, "-o",  "tpxout",ffWRITE }
 258   };
 259 #define NFILE asize(fnm)
 260
 261   /* Command line options */
 262   static real max_t = -1.0;
 263   static bool bSel=FALSE;
 264   static t_pargs pa[] = {
 265     { "-time", FALSE, etREAL, {&max_t},
 266       "Continue from frame at this time instead of the last frame" },
 267   };
 268   int nerror = 0;
 269
 270   CopyRight(stdout,argv[0]);
 271
 272   /* Parse the command line */
 273   parse_common_args(&argc,argv,0,FALSE,NFILE,fnm,asize(pa),pa,
 274                     asize(desc),desc,0,NULL);
 275
 276   bSel  = (bSel || ftp2bSet(efNDX,NFILE,fnm));
 277   bTime = opt2parg_bSet("-time",asize(pa),pa);
 278
 279   top_fn = ftp2fn(efTPX,NFILE,fnm);
 280   fprintf(stderr,"Reading toplogy and shit from %s\n",top_fn);
 281
 282   read_tpxheader(top_fn,&tpx);
 283   snew(x,tpx.natoms);
 284   snew(v,tpx.natoms);
 285   snew(ir,1);
 286   read_tpx(top_fn,&step,&run_t,&run_lambda,ir,box,&natoms,x,v,NULL,&top);
 287   run_step   = 0;
 288
 289   if (ftp2bSet(efTRN,NFILE,fnm)) {
 290     frame_fn = ftp2fn(efTRN,NFILE,fnm);
 291     fprintf(stderr,
 292             "\nREADING COORDS, VELS AND BOX FROM TRAJECTORY %s...\n\n",
 293             frame_fn);
 294
 295     fp=open_trn(frame_fn,"r");
 296     fread_trnheader(fp,&head,&bOK);
 297     if (top.atoms.nr != head.natoms)
 298       fatal_error(0,"Number of atoms in Topology (%d) "
 299                   "is not the same as in Trajectory (%d)\n",
 300                   top.atoms.nr,head.natoms);
 301     snew(newx,head.natoms);
 302     snew(newv,head.natoms);
 303     run_t      = head.t;
 304     run_step   = head.step;
 305     run_lambda = head.lambda;
 306     bOK        = fread_htrn(fp,&head,box,x,v,NULL);
 307
 308     /* Now scan until the last set of x and v (step == 0)
 309      * or the ones at step step.
 310      */
 311     bFrame = bOK;
 312     frame  = 0;
 313     while (bFrame) {
 314       fprintf(stderr,"\rRead frame %6d: step %6d time %8.3f",
 315               frame,run_step,run_t);
 316       bFrame=fread_trnheader(fp,&head,&bOK);
 317       if (bFrame || !bOK)
 318         frame++;
 319       bFrame=bFrame && bOK;
 320       if (bFrame)
 321         bOK=fread_htrn(fp,&head,newbox,newx,newv,NULL);
 322       bFrame=bFrame && bOK;
 323       if (bFrame && (head.x_size) && (head.v_size)) {
 324         tmpx = newx;
 325         newx = x;
 326         x    = tmpx;
 327         tmpv = newv;
 328         newv = v;
 329         v    = tmpv;
 330         run_t      = head.t;
 331         run_step   = head.step;
 332         run_lambda = head.lambda;
 333         copy_mat(newbox,box);
 334       }
 335       if (bTime && (head.t >= max_t))
 336         bFrame=FALSE;
 337     }
 338     close_trn(fp);
 339     fprintf(stderr,"\n");
 340     if (!bOK)
 341       fprintf(stderr,"Frame %d (step %d, time %g) is incomplete\n",
 342               frame,head.step,head.t);
 343     fprintf(stderr,"\nUsing frame of step %d time %g\n",run_step,run_t);
 344   }
 345   else {
 346     frame_fn = ftp2fn(efTPX,NFILE,fnm);
 347     fprintf(stderr,"\nUSING COORDS, VELS AND BOX FROM TPX FILE %s...\n\n",
 348             ftp2fn(efTPX,NFILE,fnm));
 349   }
 350
 351   ir->nsteps     -= run_step;
 352   ir->init_t      = run_t;
 353   ir->init_lambda = run_lambda;
 354
 355   if (!ftp2bSet(efTRN,NFILE,fnm)) {
 356     get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),1,
 357               &gnx,&index,&grpname);
 358     bSel=(gnx!=natoms);
 359     for (i=0; ((i<gnx) && (!bSel)); i++)
 360       bSel = (i!=index[i]);
 361     if (bSel) {
 362       fprintf(stderr,"Will write subset %s of original tpx containg %d "
 363               "atoms\n",grpname,gnx);
 364       reduce_topology_x(gnx,index,&top,x,v);
 365       natoms = gnx;
 366     }
 367     else
 368       fprintf(stderr,"Will write full tpx file (no selection)\n");
 369   }
 370
 371   fprintf(stderr,"Writing statusfile with starting time %g and %d steps...\n",
 372           ir->init_t,ir->nsteps);
 373   write_tpx(opt2fn("-o",NFILE,fnm),
 374             0,ir->init_t,ir->init_lambda,ir,box,
 375             natoms,x,v,NULL,&top);
 376   thanx(stdout);
 377
 378   return 0;
 379 }