src/tools/gmx_trjconv.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  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <string.h>
  40 #include <math.h>
  41 #include "macros.h"
  42 #include "assert.h"
  43 #include "sysstuff.h"
  44 #include "smalloc.h"
  45 #include "typedefs.h"
  46 #include "copyrite.h"
  47 #include "gmxfio.h"
  48 #include "tpxio.h"
  49 #include "trnio.h"
  50 #include "statutil.h"
  51 #include "futil.h"
  52 #include "pdbio.h"
  53 #include "confio.h"
  54 #include "names.h"
  55 #include "index.h"
  56 #include "vec.h"
  57 #include "xtcio.h"
  58 #include "do_fit.h"
  59 #include "rmpbc.h"
  60 #include "wgms.h"
  61 #include "magic.h"
  62 #include "pbc.h"
  63 #include "viewit.h"
  64 #include "xvgr.h"
  65 #include "gmx_ana.h"
  66
  67 enum { euSel,euRect, euTric, euCompact, euNR};
  68
  69 static real calc_isquared(int nmol,rvec m_com[],rvec m_shift[],rvec clust_com)
  70 {
  71     real Isq=0;
  72     int  i;
  73     rvec m0,dx;
  74
  75     for(i=0; (i<nmol); i++) {
  76         rvec_add(m_com[i],m_shift[i],m0);
  77         rvec_sub(m0,clust_com,dx);
  78         Isq += iprod(dx,dx);
  79     }
  80     return Isq;
  81 }
  82
  83 static void calc_pbc_cluster(int ecenter,int nrefat,t_topology *top,int ePBC,
  84                              rvec x[],atom_id index[],
  85                              rvec clust_com,matrix box)
  86 {
  87     const   real tol=1e-3;
  88     gmx_bool    bChanged;
  89     int     m,i,j,j0,j1,jj,ai,iter,is;
  90     real    fac,Isq,min_dist2;
  91     rvec    dx,ddx,xtest,xrm,box_center;
  92     int     nmol,nmol_cl,imol_center;
  93     atom_id *molind;
  94     gmx_bool    *bMol,*bTmp;
  95     rvec    *m_com,*m_shift,m0;
  96     t_pbc   pbc;
  97
  98     calc_box_center(ecenter,box,box_center);
  99
 100     /* Initiate the pbc structure */
 101     memset(&pbc,0,sizeof(pbc));
 102     set_pbc(&pbc,ePBC,box);
 103
 104     /* Convert atom index to molecular */
 105     nmol   = top->mols.nr;
 106     molind = top->mols.index;
 107     snew(bMol,nmol);
 108     snew(m_com,nmol);
 109     snew(m_shift,nmol);
 110     snew(bTmp,top->atoms.nr);
 111     nmol_cl = 0;
 112     for(i=0; (i<nrefat); i++) {
 113         /* Mark all molecules in the index */
 114         ai = index[i];
 115         bTmp[ai] = TRUE;
 116         /* Binary search assuming the molecules are sorted */
 117         j0=0;
 118         j1=nmol-1;
 119         while (j0 < j1) {
 120             if (ai < molind[j0+1])
 121                 j1 = j0;
 122             else if (ai >= molind[j1])
 123                 j0 = j1;
 124             else {
 125                 jj = (j0+j1)/2;
 126                 if (ai < molind[jj+1])
 127                     j1 = jj;
 128                 else
 129                     j0 = jj;
 130             }
 131         }
 132         bMol[j0] = TRUE;
 133     }
 134     /* Double check whether all atoms in all molecules that are marked are part
 135      * of the cluster. Simultaneously compute the center of geometry.
 136      */
 137     min_dist2   = 10*sqr(trace(box));
 138     imol_center = -1;
 139     for(i=0; (i<nmol); i++) {
 140         for(j=molind[i]; (j<molind[i+1]); j++) {
 141             if (bMol[i] && !bTmp[j])
 142                 gmx_fatal(FARGS,"Molecule %d marked for clustering but not atom %d",
 143                           i+1,j+1);
 144             else if (!bMol[i] && bTmp[j])
 145                 gmx_fatal(FARGS,"Atom %d marked for clustering but not molecule %d",
 146                           j+1,i+1);
 147             else if (bMol[i]) {
 148                 /* Compute center of geometry of molecule */
 149                 rvec_inc(m_com[i],x[j]);
 150             }
 151         }
 152         if (bMol[i]) {
 153             /* Normalize center of geometry */
 154             fac = 1.0/(molind[i+1]-molind[i]);
 155             for(m=0; (m<DIM); m++)
 156                 m_com[i][m] *= fac;
 157             /* Determine which molecule is closest to the center of the box */
 158             pbc_dx(&pbc,box_center,m_com[i],dx);
 159             if (iprod(dx,dx) < min_dist2) {
 160                 min_dist2   = iprod(dx,dx);
 161                 imol_center = i;
 162             }
 163             nmol_cl++;
 164         }
 165     }
 166     sfree(bTmp);
 167
 168     if (nmol_cl <= 0) {
 169         fprintf(stderr,"No molecules selected in the cluster\n");
 170         return;
 171     } else if (imol_center == -1) {
 172         fprintf(stderr,"No central molecules could be found\n");
 173         return;
 174     }
 175
 176
 177     /* First calculation is incremental */
 178     clear_rvec(clust_com);
 179     Isq = 0;
 180     for(i=m=0; (i<nmol); i++) {
 181         /* Check whether this molecule is part of the cluster */
 182         if (bMol[i]) {
 183             if ((i > 0) && (m > 0)) {
 184                 /* Compute center of cluster by dividing by number of molecules */
 185                 svmul(1.0/m,clust_com,xrm);
 186                 /* Distance vector between molecular COM and cluster */
 187                 pbc_dx(&pbc,m_com[i],xrm,dx);
 188                 rvec_add(xrm,dx,xtest);
 189                 /* xtest is now the image of m_com[i] that is closest to clust_com */
 190                 rvec_inc(clust_com,xtest);
 191                 rvec_sub(xtest,m_com[i],m_shift[i]);
 192             }
 193             else {
 194                 rvec_inc(clust_com,m_com[i]);
 195             }
 196             m++;
 197         }
 198     }
 199     assert(m == nmol_cl);
 200     svmul(1/nmol_cl,clust_com,clust_com);
 201     put_atom_in_box(box,clust_com);
 202
 203     /* Now check if any molecule is more than half the box from the COM */
 204     if (box) {
 205         iter = 0;
 206         do {
 207             bChanged = FALSE;
 208             for(i=0; (i<nmol) && !bChanged; i++) {
 209                 if (bMol[i]) {
 210                     /* Sum com and shift from com */
 211                     rvec_add(m_com[i],m_shift[i],m0);
 212                     pbc_dx(&pbc,m0,clust_com,dx);
 213                     rvec_add(clust_com,dx,xtest);
 214                     rvec_sub(xtest,m0,ddx);
 215                     if (iprod(ddx,ddx) > tol) {
 216                         /* Here we have used the wrong image for contributing to the COM */
 217                         rvec_sub(xtest,m_com[i],m_shift[i]);
 218                         for(j=0; (j<DIM); j++)
 219                             clust_com[j] += (xtest[j]-m0[j])/nmol_cl;
 220                         bChanged = TRUE;
 221                     }
 222                 }
 223             }
 224             Isq = calc_isquared(nmol,m_com,m_shift,clust_com);
 225             put_atom_in_box(box,clust_com);
 226
 227             if (bChanged && (iter > 0))
 228                 printf("COM: %8.3f  %8.3f  %8.3f  iter = %d  Isq = %8.3f\n",
 229                        clust_com[XX],clust_com[YY],clust_com[ZZ],iter,Isq);
 230             iter++;
 231         } while (bChanged);
 232     }
 233     /* Now transfer the shift to all atoms in the molecule */
 234     for(i=0; (i<nmol); i++) {
 235         if (bMol[i]) {
 236             for(j=molind[i]; (j<molind[i+1]); j++)
 237                 rvec_inc(x[j],m_shift[i]);
 238         }
 239     }
 240     sfree(bMol);
 241     sfree(m_com);
 242     sfree(m_shift);
 243 }
 244
 245 static void put_molecule_com_in_box(int unitcell_enum,int ecenter,
 246                                     t_block *mols,
 247                                     int natoms,t_atom atom[],
 248                                     int ePBC,matrix box,rvec x[])
 249 {
 250     atom_id i,j;
 251     int     d;
 252     rvec    com,new_com,shift,dx,box_center;
 253     real    m;
 254     double  mtot;
 255     t_pbc   pbc;
 256
 257     calc_box_center(ecenter,box,box_center);
 258     set_pbc(&pbc,ePBC,box);
 259     if (mols->nr <= 0)
 260         gmx_fatal(FARGS,"There are no molecule descriptions. I need a tpr file for this pbc option.");
 261     for(i=0; (i<mols->nr); i++) {
 262         /* calc COM */
 263         clear_rvec(com);
 264         mtot = 0;
 265         for(j=mols->index[i]; (j<mols->index[i+1] && j<natoms); j++) {
 266             m = atom[j].m;
 267             for(d=0; d<DIM; d++)
 268                 com[d] += m*x[j][d];
 269             mtot += m;
 270         }
 271         /* calculate final COM */
 272         svmul(1.0/mtot, com, com);
 273
 274         /* check if COM is outside box */
 275         copy_rvec(com,new_com);
 276         switch (unitcell_enum) {
 277         case euRect:
 278             put_atoms_in_box(box,1,&new_com);
 279             break;
 280         case euTric:
 281             put_atoms_in_triclinic_unitcell(ecenter,box,1,&new_com);
 282             break;
 283         case euCompact:
 284             put_atoms_in_compact_unitcell(ePBC,ecenter,box,1,&new_com);
 285             break;
 286         }
 287         rvec_sub(new_com,com,shift);
 288         if (norm2(shift) > 0) {
 289             if (debug)
 290                 fprintf (debug,"\nShifting position of molecule %d "
 291                          "by %8.3f  %8.3f  %8.3f\n", i+1, PR_VEC(shift));
 292             for(j=mols->index[i]; (j<mols->index[i+1] && j<natoms); j++) {
 293                 rvec_inc(x[j],shift);
 294             }
 295         }
 296     }
 297 }
 298
 299 static void put_residue_com_in_box(int unitcell_enum,int ecenter,
 300                                    int natoms,t_atom atom[],
 301                                    int ePBC,matrix box,rvec x[])
 302 {
 303     atom_id i, j, res_start, res_end, res_nat;
 304     int     d, presnr;
 305     real    m;
 306     double  mtot;
 307     rvec    box_center,com,new_com,shift;
 308
 309     calc_box_center(ecenter,box,box_center);
 310
 311     presnr = NOTSET;
 312     res_start = 0;
 313     clear_rvec(com);
 314     mtot = 0;
 315     for(i=0; i<natoms+1; i++) {
 316         if (i == natoms || (presnr != atom[i].resind && presnr != NOTSET)) {
 317             /* calculate final COM */
 318             res_end = i;
 319             res_nat = res_end - res_start;
 320             svmul(1.0/mtot, com, com);
 321
 322             /* check if COM is outside box */
 323             copy_rvec(com,new_com);
 324             switch (unitcell_enum) {
 325             case euRect:
 326                 put_atoms_in_box(box,1,&new_com);
 327                 break;
 328             case euTric:
 329                 put_atoms_in_triclinic_unitcell(ecenter,box,1,&new_com);
 330                 break;
 331             case euCompact:
 332                 put_atoms_in_compact_unitcell(ePBC,ecenter,box,1,&new_com);
 333                 break;
 334             }
 335             rvec_sub(new_com,com,shift);
 336             if (norm2(shift)) {
 337                 if (debug)
 338                     fprintf (debug,"\nShifting position of residue %d (atoms %u-%u) "
 339                              "by %g,%g,%g\n", atom[res_start].resind+1,
 340                              res_start+1, res_end+1, PR_VEC(shift));
 341                 for(j=res_start; j<res_end; j++)
 342                     rvec_inc(x[j],shift);
 343             }
 344             clear_rvec(com);
 345             mtot = 0;
 346
 347             /* remember start of new residue */
 348             res_start = i;
 349         }
 350         if (i < natoms) {
 351             /* calc COM */
 352             m = atom[i].m;
 353             for(d=0; d<DIM; d++)
 354                 com[d] += m*x[i][d];
 355             mtot += m;
 356
 357             presnr = atom[i].resind;
 358         }
 359     }
 360 }
 361
 362 static void center_x(int ecenter,rvec x[],matrix box,int n,int nc,atom_id ci[])
 363 {
 364     int i,m,ai;
 365     rvec cmin,cmax,box_center,dx;
 366
 367     if (nc > 0) {
 368         copy_rvec(x[ci[0]],cmin);
 369         copy_rvec(x[ci[0]],cmax);
 370         for(i=0; i<nc; i++) {
 371             ai=ci[i];
 372             for(m=0; m<DIM; m++) {
 373                 if (x[ai][m] < cmin[m])
 374                     cmin[m]=x[ai][m];
 375                 else if (x[ai][m] > cmax[m])
 376                     cmax[m]=x[ai][m];
 377             }
 378         }
 379         calc_box_center(ecenter,box,box_center);
 380         for(m=0; m<DIM; m++)
 381             dx[m] = box_center[m]-(cmin[m]+cmax[m])*0.5;
 382
 383         for(i=0; i<n; i++)
 384             rvec_inc(x[i],dx);
 385     }
 386 }
 387
 388 static void mk_filenm(char *base,const char *ext,int ndigit,int file_nr,
 389                       char out_file[])
 390 {
 391     char nbuf[128];
 392     int  nd=0,fnr;
 393
 394     strcpy(out_file,base);
 395     fnr = file_nr;
 396     do {
 397         fnr /= 10;
 398         nd++;
 399     } while (fnr > 0);
 400
 401     if (nd < ndigit)
 402         strncat(out_file,"00000000000",ndigit-nd);
 403     sprintf(nbuf,"%d.",file_nr);
 404     strcat(out_file,nbuf);
 405     strcat(out_file,ext);
 406 }
 407
 408 void check_trn(const char *fn)
 409 {
 410     if ((fn2ftp(fn) != efTRJ)  && (fn2ftp(fn) != efTRR))
 411         gmx_fatal(FARGS,"%s is not a trj file, exiting\n",fn);
 412 }
 413
 414 #if (!defined WIN32 && !defined _WIN32 && !defined WIN64 && !defined _WIN64)
 415 void do_trunc(const char *fn, real t0)
 416 {
 417     t_fileio     *in;
 418     FILE         *fp;
 419     gmx_bool         bStop,bOK;
 420     t_trnheader  sh;
 421     gmx_off_t    fpos;
 422     char         yesno[256];
 423     int          j;
 424     real         t=0;
 425
 426     if (t0 == -1)
 427         gmx_fatal(FARGS,"You forgot to set the truncation time");
 428
 429     /* Check whether this is a .trj file */
 430     check_trn(fn);
 431
 432     in   = open_trn(fn,"r");
 433     fp   = gmx_fio_getfp(in);
 434     if (fp == NULL) {
 435         fprintf(stderr,"Sorry, can not trunc %s, truncation of this filetype is not supported\n",fn);
 436         close_trn(in);
 437     } else {
 438         j    = 0;
 439         fpos = gmx_fio_ftell(in);
 440         bStop= FALSE;
 441         while (!bStop && fread_trnheader(in,&sh,&bOK)) {
 442             fread_htrn(in,&sh,NULL,NULL,NULL,NULL);
 443             fpos=gmx_ftell(fp);
 444             t=sh.t;
 445             if (t>=t0) {
 446                 gmx_fseek(fp, fpos, SEEK_SET);
 447                 bStop=TRUE;
 448             }
 449         }
 450         if (bStop) {
 451             fprintf(stderr,"Do you REALLY want to truncate this trajectory (%s) at:\n"
 452                     "frame %d, time %g, bytes %ld ??? (type YES if so)\n",
 453                     fn,j,t,(long int)fpos);
 454             if(1 != scanf("%s",yesno))
 455             {
 456                 gmx_fatal(FARGS,"Error reading user input");
 457             }
 458             if (strcmp(yesno,"YES") == 0) {
 459                 fprintf(stderr,"Once again, I'm gonna DO this...\n");
 460                 close_trn(in);
 461                 if(0 != truncate(fn,fpos))
 462                 {
 463                     gmx_fatal(FARGS,"Error truncating file %s",fn);
 464                 }
 465             }
 466             else {
 467                 fprintf(stderr,"Ok, I'll forget about it\n");
 468             }
 469         }
 470         else {
 471             fprintf(stderr,"Already at end of file (t=%g)...\n",t);
 472             close_trn(in);
 473         }
 474     }
 475 }
 476 #endif
 477
 478 int gmx_trjconv(int argc,char *argv[])
 479 {
 480     const char *desc[] = {
 481         "trjconv can convert trajectory files in many ways:[BR]",
 482         "[BB]1.[bb] from one format to another[BR]",
 483         "[BB]2.[bb] select a subset of atoms[BR]",
 484         "[BB]3.[bb] change the periodicity representation[BR]",
 485         "[BB]4.[bb] keep multimeric molecules together[BR]",
 486         "[BB]5.[bb] center atoms in the box[BR]",
 487         "[BB]6.[bb] fit atoms to reference structure[BR]",
 488         "[BB]7.[bb] reduce the number of frames[BR]",
 489         "[BB]8.[bb] change the timestamps of the frames ",
 490         "([TT]-t0[tt] and [TT]-timestep[tt])[BR]",
 491         "[BB]9.[bb] cut the trajectory in small subtrajectories according",
 492         "to information in an index file. This allows subsequent analysis of",
 493         "the subtrajectories that could, for example be the result of a",
 494         "cluster analysis. Use option [TT]-sub[tt].",
 495         "This assumes that the entries in the index file are frame numbers and",
 496         "dumps each group in the index file to a separate trajectory file.[BR]",
 497         "[BB]10.[bb] select frames within a certain range of a quantity given",
 498         "in an [TT].xvg[tt] file.[PAR]",
 499         "The program [TT]trjcat[tt] can concatenate multiple trajectory files.",
 500         "[PAR]",
 501         "Currently seven formats are supported for input and output:",
 502         "[TT].xtc[tt], [TT].trr[tt], [TT].trj[tt], [TT].gro[tt], [TT].g96[tt],",
 503         "[TT].pdb[tt] and [TT].g87[tt].",
 504         "The file formats are detected from the file extension.",
 505         "The precision of [TT].xtc[tt] and [TT].gro[tt] output is taken from the",
 506         "input file for [TT].xtc[tt], [TT].gro[tt] and [TT].pdb[tt],",
 507         "and from the [TT]-ndec[tt] option for other input formats. The precision",
 508         "is always taken from [TT]-ndec[tt], when this option is set.",
 509         "All other formats have fixed precision. [TT].trr[tt] and [TT].trj[tt]",
 510         "output can be single or double precision, depending on the precision",
 511         "of the trjconv binary.",
 512         "Note that velocities are only supported in",
 513         "[TT].trr[tt], [TT].trj[tt], [TT].gro[tt] and [TT].g96[tt] files.[PAR]",
 514         "Option [TT]-app[tt] can be used to",
 515         "append output to an existing trajectory file.",
 516         "No checks are performed to ensure integrity",
 517         "of the resulting combined trajectory file.[PAR]",
 518         "Option [TT]-sep[tt] can be used to write every frame to a separate",
 519         ".gro, .g96 or .pdb file, default all frames all written to one file.",
 520         "[TT].pdb[tt] files with all frames concatenated can be viewed with",
 521         "[TT]rasmol -nmrpdb[tt].[PAR]",
 522         "It is possible to select part of your trajectory and write it out",
 523         "to a new trajectory file in order to save disk space, e.g. for leaving",
 524         "out the water from a trajectory of a protein in water.",
 525         "[BB]ALWAYS[bb] put the original trajectory on tape!",
 526         "We recommend to use the portable [TT].xtc[tt] format for your analysis",
 527         "to save disk space and to have portable files.[PAR]",
 528         "There are two options for fitting the trajectory to a reference",
 529         "either for essential dynamics analysis or for whatever.",
 530         "The first option is just plain fitting to a reference structure",
 531         "in the structure file, the second option is a progressive fit",
 532         "in which the first timeframe is fitted to the reference structure ",
 533         "in the structure file to obtain and each subsequent timeframe is ",
 534         "fitted to the previously fitted structure. This way a continuous",
 535         "trajectory is generated, which might not be the case when using the",
 536         "regular fit method, e.g. when your protein undergoes large",
 537         "conformational transitions.[PAR]",
 538         "Option [TT]-pbc[tt] sets the type of periodic boundary condition",
 539         "treatment:[BR]",
 540         "* [TT]mol[tt] puts the center of mass of molecules in the box.[BR]",
 541         "* [TT]res[tt] puts the center of mass of residues in the box.[BR]",
 542         "* [TT]atom[tt] puts all the atoms in the box.[BR]",
 543         "* [TT]nojump[tt] checks if atoms jump across the box and then puts",
 544         "them back. This has the effect that all molecules",
 545         "will remain whole (provided they were whole in the initial",
 546         "conformation), note that this ensures a continuous trajectory but",
 547         "molecules may diffuse out of the box. The starting configuration",
 548         "for this procedure is taken from the structure file, if one is",
 549         "supplied, otherwise it is the first frame.[BR]",
 550         "* [TT]cluster[tt] clusters all the atoms in the selected index",
 551         "such that they are all closest to the center of mass of the cluster",
 552         "which is iteratively updated. Note that this will only give meaningful",
 553         "results if you in fact have a cluster. Luckily that can be checked",
 554         "afterwards using a trajectory viewer. Note also that if your molecules",
 555         "are broken this will not work either.[BR]",
 556         "* [TT]whole[tt] only makes broken molecules whole.[PAR]",
 557         "Option [TT]-ur[tt] sets the unit cell representation for options",
 558         "[TT]mol[tt], [TT]res[tt] and [TT]atom[tt] of [TT]-pbc[tt].",
 559         "All three options give different results for triclinic boxes and",
 560         "identical results for rectangular boxes.",
 561         "[TT]rect[tt] is the ordinary brick shape.",
 562         "[TT]tric[tt] is the triclinic unit cell.",
 563         "[TT]compact[tt] puts all atoms at the closest distance from the center",
 564         "of the box. This can be useful for visualizing e.g. truncated",
 565         "octahedrons. The center for options [TT]tric[tt] and [TT]compact[tt]",
 566         "is [TT]tric[tt] (see below), unless the option [TT]-boxcenter[tt]",
 567         "is set differently.[PAR]",
 568         "Option [TT]-center[tt] centers the system in the box. The user can",
 569         "select the group which is used to determine the geometrical center.",
 570         "Option [TT]-boxcenter[tt] sets the location of the center of the box",
 571         "for options [TT]-pbc[tt] and [TT]-center[tt]. The center options are:",
 572         "[TT]tric[tt]: half of the sum of the box vectors,",
 573         "[TT]rect[tt]: half of the box diagonal,",
 574         "[TT]zero[tt]: zero.",
 575         "Use option [TT]-pbc mol[tt] in addition to [TT]-center[tt] when you",
 576         "want all molecules in the box after the centering.[PAR]",
 577         "With [TT]-dt[tt] it is possible to reduce the number of ",
 578         "frames in the output. This option relies on the accuracy of the times",
 579         "in your input trajectory, so if these are inaccurate use the",
 580         "[TT]-timestep[tt] option to modify the time (this can be done",
 581         "simultaneously). For making smooth movies the program [TT]g_filter[tt]",
 582         "can reduce the number of frames while using low-pass frequency",
 583         "filtering, this reduces aliasing of high frequency motions.[PAR]",
 584         "Using [TT]-trunc[tt] trjconv can truncate [TT].trj[tt] in place, i.e.",
 585         "without copying the file. This is useful when a run has crashed",
 586         "during disk I/O (one more disk full), or when two contiguous",
 587         "trajectories must be concatenated without have double frames.[PAR]",
 588         "[TT]trjcat[tt] is more suitable for concatenating trajectory files.[PAR]",
 589         "Option [TT]-dump[tt] can be used to extract a frame at or near",
 590         "one specific time from your trajectory.[PAR]",
 591         "Option [TT]-drop[tt] reads an [TT].xvg[tt] file with times and values.",
 592         "When options [TT]-dropunder[tt] and/or [TT]-dropover[tt] are set,",
 593         "frames with a value below and above the value of the respective options",
 594         "will not be written."
 595     };
 596
 597     int pbc_enum;
 598     enum
 599     {
 600         epSel,
 601         epNone,
 602         epComMol,
 603         epComRes,
 604         epComAtom,
 605         epNojump,
 606         epCluster,
 607         epWhole,
 608         epNR
 609     };
 610     const char *pbc_opt[epNR + 1] =
 611         { NULL, "none", "mol", "res", "atom", "nojump", "cluster", "whole",
 612             NULL };
 613
 614     int unitcell_enum;
 615     const char *unitcell_opt[euNR+1] =
 616         { NULL, "rect", "tric", "compact", NULL };
 617
 618     enum
 619     { ecSel, ecTric, ecRect, ecZero, ecNR};
 620     const char *center_opt[ecNR+1] =
 621         { NULL, "tric", "rect", "zero", NULL };
 622     int ecenter;
 623
 624     int fit_enum;
 625     enum
 626     {
 627         efSel, efNone, efFit, efFitXY, efReset, efResetXY, efPFit, efNR
 628     };
 629     const char *fit[efNR + 1] =
 630         { NULL, "none", "rot+trans", "rotxy+transxy", "translation", "transxy",
 631             "progressive", NULL };
 632
 633     static gmx_bool  bAppend=FALSE,bSeparate=FALSE,bVels=TRUE,bForce=FALSE,bCONECT=FALSE;
 634     static gmx_bool  bCenter=FALSE;
 635     static int   skip_nr=1,ndec=3,nzero=0;
 636     static real  tzero=0,delta_t=0,timestep=0,ttrunc=-1,tdump=-1,split_t=0;
 637     static rvec  newbox = {0,0,0}, shift = {0,0,0}, trans = {0,0,0};
 638     static char  *exec_command=NULL;
 639     static real  dropunder=0,dropover=0;
 640     static gmx_bool  bRound=FALSE;
 641
 642     t_pargs
 643         pa[] =
 644             {
 645                     { "-skip", FALSE, etINT,
 646                         { &skip_nr }, "Only write every nr-th frame" },
 647                     { "-dt", FALSE, etTIME,
 648                         { &delta_t },
 649                         "Only write frame when t MOD dt = first time (%t)" },
 650                     { "-round", FALSE, etBOOL,
 651                         { &bRound }, "Round measurements to nearest picosecond"
 652                     },
 653                     { "-dump", FALSE, etTIME,
 654                         { &tdump }, "Dump frame nearest specified time (%t)" },
 655                     { "-t0", FALSE, etTIME,
 656                         { &tzero },
 657                         "Starting time (%t) (default: don't change)" },
 658                     { "-timestep", FALSE, etTIME,
 659                         { &timestep },
 660                         "Change time step between input frames (%t)" },
 661                     { "-pbc", FALSE, etENUM,
 662                         { pbc_opt },
 663                         "PBC treatment (see help text for full description)" },
 664                     { "-ur", FALSE, etENUM,
 665                         { unitcell_opt }, "Unit-cell representation" },
 666                     { "-center", FALSE, etBOOL,
 667                         { &bCenter }, "Center atoms in box" },
 668                     { "-boxcenter", FALSE, etENUM,
 669                         { center_opt }, "Center for -pbc and -center" },
 670                     { "-box", FALSE, etRVEC,
 671                         { newbox },
 672                         "Size for new cubic box (default: read from input)" },
 673                     { "-trans", FALSE, etRVEC,
 674                         { trans },
 675                         "All coordinates will be translated by trans. This "
 676                         "can advantageously be combined with -pbc mol -ur "
 677                         "compact." },
 678                     { "-shift", FALSE, etRVEC,
 679                         { shift },
 680                         "All coordinates will be shifted by framenr*shift" },
 681                     { "-fit", FALSE, etENUM,
 682                         { fit },
 683                         "Fit molecule to ref structure in the structure file" },
 684                     { "-ndec", FALSE, etINT,
 685                         { &ndec },
 686                         "Precision for .xtc and .gro writing in number of "
 687                         "decimal places" },
 688                     { "-vel", FALSE, etBOOL,
 689                         { &bVels }, "Read and write velocities if possible" },
 690                     { "-force", FALSE, etBOOL,
 691                         { &bForce }, "Read and write forces if possible" },
 692 #if (!defined WIN32 && !defined _WIN32 && !defined WIN64 && !defined _WIN64)
 693                     { "-trunc", FALSE, etTIME,
 694                         { &ttrunc },
 695                         "Truncate input trj file after this time (%t)" },
 696 #endif
 697                     { "-exec", FALSE, etSTR,
 698                         { &exec_command },
 699                         "Execute command for every output frame with the "
 700                         "frame number as argument" },
 701                     { "-app", FALSE, etBOOL,
 702                         { &bAppend }, "Append output" },
 703                     { "-split", FALSE, etTIME,
 704                         { &split_t },
 705                         "Start writing new file when t MOD split = first "
 706                         "time (%t)" },
 707                     { "-sep", FALSE, etBOOL,
 708                         { &bSeparate },
 709                         "Write each frame to a separate .gro, .g96 or .pdb "
 710                         "file" },
 711                     { "-nzero", FALSE, etINT,
 712                         { &nzero },
 713                         "If the -sep flag is set, use these many digits "
 714                         "for the file numbers and prepend zeros as needed" },
 715                     { "-dropunder", FALSE, etREAL,
 716                         { &dropunder }, "Drop all frames below this value" },
 717                     { "-dropover", FALSE, etREAL,
 718                         { &dropover }, "Drop all frames above this value" },
 719                     { "-conect", FALSE, etBOOL,
 720                         { &bCONECT },
 721                         "Add conect records when writing pdb files. Useful "
 722                         "for visualization of non-standard molecules, e.g. "
 723                         "coarse grained ones" } };
 724 #define NPA asize(pa)
 725
 726     FILE         *out=NULL;
 727     t_trxstatus *trxout=NULL;
 728     t_trxstatus *status;
 729     int          ftp,ftpin=0,file_nr;
 730     t_trxframe   fr,frout;
 731     int          flags;
 732     rvec         *xmem=NULL,*vmem=NULL,*fmem=NULL;
 733     rvec         *xp=NULL,x_shift,hbox,box_center,dx;
 734     real         xtcpr, lambda,*w_rls=NULL;
 735     int          m,i,d,frame,outframe,natoms,nout,ncent,nre,newstep=0,model_nr;
 736 #define SKIP 10
 737     t_topology   top;
 738     gmx_conect   gc=NULL;
 739     int          ePBC=-1;
 740     t_atoms      *atoms=NULL,useatoms;
 741     matrix       top_box;
 742     atom_id      *index,*cindex;
 743     char         *grpnm;
 744     int          *frindex,nrfri;
 745     char         *frname;
 746     int          ifit,irms,my_clust=-1;
 747     atom_id      *ind_fit,*ind_rms;
 748     char         *gn_fit,*gn_rms;
 749     t_cluster_ndx *clust=NULL;
 750     t_trxstatus  **clust_status=NULL;
 751     int          *clust_status_id=NULL;
 752     int          ntrxopen=0;
 753     int          *nfwritten=NULL;
 754     int          ndrop=0,ncol,drop0=0,drop1=0,dropuse=0;
 755     double       **dropval;
 756     real         tshift=0,t0=-1,dt=0.001,prec;
 757     gmx_bool         bFit,bFitXY,bPFit,bReset;
 758     int          nfitdim;
 759     gmx_rmpbc_t  gpbc=NULL;
 760     gmx_bool         bRmPBC,bPBCWhole,bPBCcomRes,bPBCcomMol,bPBCcomAtom,bPBC,bNoJump,bCluster;
 761     gmx_bool         bCopy,bDoIt,bIndex,bTDump,bSetTime,bTPS=FALSE,bDTset=FALSE;
 762     gmx_bool         bExec,bTimeStep=FALSE,bDumpFrame=FALSE,bSetPrec,bNeedPrec;
 763     gmx_bool         bHaveFirstFrame,bHaveNextFrame,bSetBox,bSetUR,bSplit=FALSE;
 764     gmx_bool         bSubTraj=FALSE,bDropUnder=FALSE,bDropOver=FALSE,bTrans=FALSE;
 765     gmx_bool         bWriteFrame,bSplitHere;
 766     const char   *top_file,*in_file,*out_file=NULL;
 767     char         out_file2[256],*charpt;
 768     char         *outf_base=NULL;
 769     const char   *outf_ext=NULL;
 770     char         top_title[256],title[256],command[256],filemode[5];
 771     int          xdr=0;
 772     gmx_bool         bWarnCompact=FALSE;
 773     const char  *warn;
 774     output_env_t oenv;
 775
 776     t_filenm fnm[] = {
 777         { efTRX, "-f",   NULL,      ffREAD  },
 778         { efTRO, "-o",   NULL,      ffWRITE },
 779         { efTPS, NULL,   NULL,      ffOPTRD },
 780         { efNDX, NULL,   NULL,      ffOPTRD },
 781         { efNDX, "-fr",  "frames",  ffOPTRD },
 782         { efNDX, "-sub", "cluster", ffOPTRD },
 783         { efXVG, "-drop","drop",    ffOPTRD }
 784     };
 785 #define NFILE asize(fnm)
 786
 787     CopyRight(stderr,argv[0]);
 788     parse_common_args(&argc,argv,
 789                       PCA_CAN_BEGIN | PCA_CAN_END | PCA_CAN_VIEW |
 790                       PCA_TIME_UNIT | PCA_BE_NICE,
 791                       NFILE,fnm,NPA,pa,asize(desc),desc,
 792                       0,NULL,&oenv);
 793
 794     top_file = ftp2fn(efTPS,NFILE,fnm);
 795     init_top(&top);
 796
 797     /* Check command line */
 798     in_file=opt2fn("-f",NFILE,fnm);
 799
 800     if (ttrunc != -1) {
 801 #if (!defined WIN32 && !defined _WIN32 && !defined WIN64 && !defined _WIN64)
 802         do_trunc(in_file,ttrunc);
 803 #endif
 804     }
 805     else {
 806         /* mark active cmdline options */
 807         bSetBox   = opt2parg_bSet("-box", NPA, pa);
 808         bSetTime  = opt2parg_bSet("-t0", NPA, pa);
 809         bSetPrec  = opt2parg_bSet("-ndec", NPA, pa);
 810         bSetUR    = opt2parg_bSet("-ur", NPA, pa);
 811         bExec     = opt2parg_bSet("-exec", NPA, pa);
 812         bTimeStep = opt2parg_bSet("-timestep", NPA, pa);
 813         bTDump    = opt2parg_bSet("-dump", NPA, pa);
 814         bDropUnder= opt2parg_bSet("-dropunder", NPA, pa);
 815         bDropOver = opt2parg_bSet("-dropover", NPA, pa);
 816         bTrans    = opt2parg_bSet("-trans",NPA,pa);
 817         bSplit    = (split_t != 0);
 818
 819         /* parse enum options */
 820         fit_enum   = nenum(fit);
 821         bFit       = (fit_enum==efFit || fit_enum==efFitXY);
 822         bFitXY     = fit_enum==efFitXY;
 823         bReset     = (fit_enum==efReset || fit_enum==efResetXY);
 824         bPFit      = fit_enum==efPFit;
 825         pbc_enum   = nenum(pbc_opt);
 826         bPBCWhole  = pbc_enum==epWhole;
 827         bPBCcomRes = pbc_enum==epComRes;
 828         bPBCcomMol = pbc_enum==epComMol;
 829         bPBCcomAtom= pbc_enum==epComAtom;
 830         bNoJump    = pbc_enum==epNojump;
 831         bCluster   = pbc_enum==epCluster;
 832         bPBC       = pbc_enum!=epNone;
 833         unitcell_enum = nenum(unitcell_opt);
 834         ecenter    = nenum(center_opt) - ecTric;
 835
 836         /* set and check option dependencies */
 837         if (bPFit) bFit = TRUE; /* for pfit, fit *must* be set */
 838         if (bFit) bReset = TRUE; /* for fit, reset *must* be set */
 839         nfitdim = 0;
 840         if (bFit || bReset)
 841             nfitdim = (fit_enum==efFitXY || fit_enum==efResetXY) ? 2 : 3;
 842         bRmPBC = bFit || bPBCWhole || bPBCcomRes || bPBCcomMol;
 843
 844         if (bSetUR) {
 845             if (!(bPBCcomRes || bPBCcomMol ||  bPBCcomAtom)) {
 846                 fprintf(stderr,
 847                         "WARNING: Option for unitcell representation (-ur %s)\n"
 848                         "         only has effect in combination with -pbc %s, %s or %s.\n"
 849                         "         Ingoring unitcell representation.\n\n",
 850                         unitcell_opt[0],pbc_opt[2],pbc_opt[3],pbc_opt[4]);
 851                 bSetUR = FALSE;
 852             }
 853         }
 854         if (bFit && bPBC) {
 855             gmx_fatal(FARGS,"PBC condition treatment does not work together with rotational fit.\n"
 856                       "Please do the PBC condition treatment first and then run trjconv in a second step\n"
 857                       "for the rotational fit.\n"
 858                       "First doing the rotational fit and then doing the PBC treatment gives incorrect\n"
 859                       "results!");
 860         }
 861
 862         /* ndec is in nr of decimal places, prec is a multiplication factor: */
 863         prec = 1;
 864         for (i=0; i<ndec; i++)
 865             prec *= 10;
 866
 867         bIndex=ftp2bSet(efNDX,NFILE,fnm);
 868
 869
 870         /* Determine output type */
 871         out_file=opt2fn("-o",NFILE,fnm);
 872         ftp=fn2ftp(out_file);
 873         fprintf(stderr,"Will write %s: %s\n",ftp2ext(ftp),ftp2desc(ftp));
 874         bNeedPrec = (ftp==efXTC || ftp==efGRO);
 875         if (bVels) {
 876             /* check if velocities are possible in input and output files */
 877             ftpin=fn2ftp(in_file);
 878             bVels= (ftp==efTRR || ftp==efTRJ || ftp==efGRO || ftp==efG96)
 879                 && (ftpin==efTRR || ftpin==efTRJ || ftpin==efGRO || ftpin==efG96 ||
 880                     ftpin==efCPT);
 881         }
 882         if (bSeparate || bSplit) {
 883             outf_ext = strrchr(out_file,'.');
 884             if (outf_ext == NULL)
 885                 gmx_fatal(FARGS,"Output file name '%s' does not contain a '.'",out_file);
 886             outf_base = strdup(out_file);
 887             outf_base[outf_ext - out_file] = '\0';
 888         }
 889
 890         bSubTraj = opt2bSet("-sub",NFILE,fnm);
 891         if (bSubTraj) {
 892             if ((ftp != efXTC) && (ftp != efTRN))
 893                 gmx_fatal(FARGS,"Can only use the sub option with output file types "
 894                           "xtc and trr");
 895             clust = cluster_index(NULL,opt2fn("-sub",NFILE,fnm));
 896
 897             /* Check for number of files disabled, as FOPEN_MAX is not the correct
 898              * number to check for. In my linux box it is only 16.
 899              */
 900             if (0 && (clust->clust->nr > FOPEN_MAX-4))
 901                 gmx_fatal(FARGS,"Can not open enough (%d) files to write all the"
 902                           " trajectories.\ntry splitting the index file in %d parts.\n"
 903                           "FOPEN_MAX = %d",
 904                           clust->clust->nr,1+clust->clust->nr/FOPEN_MAX,FOPEN_MAX);
 905
 906             snew(clust_status,clust->clust->nr);
 907             snew(clust_status_id,clust->clust->nr);
 908             snew(nfwritten,clust->clust->nr);
 909             for(i=0; (i<clust->clust->nr); i++)
 910             {
 911                 clust_status[i] = NULL;
 912                 clust_status_id[i] = -1;
 913             }
 914             bSeparate = bSplit = FALSE;
 915         }
 916         /* skipping */
 917         if (skip_nr <= 0) {
 918         }
 919
 920         /* Determine whether to read a topology */
 921         bTPS = (ftp2bSet(efTPS,NFILE,fnm) ||
 922             bRmPBC || bReset || bPBCcomMol || bCluster ||
 923             (ftp == efGRO) || (ftp == efPDB) || bCONECT);
 924
 925         /* Determine if when can read index groups */
 926         bIndex = (bIndex || bTPS);
 927
 928         if (bTPS) {
 929             read_tps_conf(top_file,top_title,&top,&ePBC,&xp,NULL,top_box,
 930                           bReset || bPBCcomRes);
 931             atoms=&top.atoms;
 932
 933             if (0 == top.mols.nr && (bCluster || bPBCcomMol))
 934             {
 935                 gmx_fatal(FARGS,"Option -pbc %s requires a .tpr file for the -s option",pbc_opt[pbc_enum]);
 936             }
 937
 938             /* top_title is only used for gro and pdb,
 939              * the header in such a file is top_title t= ...
 940              * to prevent a double t=, remove it from top_title
 941              */
 942             if ((charpt=strstr(top_title," t= ")))
 943                 charpt[0]='\0';
 944
 945             if (bCONECT)
 946                 gc = gmx_conect_generate(&top);
 947             if (bRmPBC)
 948               gpbc = gmx_rmpbc_init(&top.idef,ePBC,top.atoms.nr,top_box);
 949         }
 950
 951         /* get frame number index */
 952         frindex=NULL;
 953         if (opt2bSet("-fr",NFILE,fnm)) {
 954             printf("Select groups of frame number indices:\n");
 955             rd_index(opt2fn("-fr",NFILE,fnm),1,&nrfri,(atom_id **)&frindex,&frname);
 956             if (debug)
 957                 for(i=0; i<nrfri; i++)
 958                     fprintf(debug,"frindex[%4d]=%4d\n",i,frindex[i]);
 959         }
 960
 961         /* get index groups etc. */
 962         if (bReset) {
 963             printf("Select group for %s fit\n",
 964                    bFit?"least squares":"translational");
 965             get_index(atoms,ftp2fn_null(efNDX,NFILE,fnm),
 966                       1,&ifit,&ind_fit,&gn_fit);
 967
 968             if (bFit) {
 969                 if (ifit < 2)
 970                     gmx_fatal(FARGS,"Need at least 2 atoms to fit!\n");
 971                 else if (ifit == 3)
 972                     fprintf(stderr,"WARNING: fitting with only 2 atoms is not unique\n");
 973             }
 974         }
 975         else if (bCluster) {
 976             printf("Select group for clustering\n");
 977             get_index(atoms,ftp2fn_null(efNDX,NFILE,fnm),
 978                       1,&ifit,&ind_fit,&gn_fit);
 979         }
 980
 981         if (bIndex) {
 982             if (bCenter) {
 983                 printf("Select group for centering\n");
 984                 get_index(atoms,ftp2fn_null(efNDX,NFILE,fnm),
 985                           1,&ncent,&cindex,&grpnm);
 986             }
 987             printf("Select group for output\n");
 988             get_index(atoms,ftp2fn_null(efNDX,NFILE,fnm),
 989                       1,&nout,&index,&grpnm);
 990         } else {
 991             /* no index file, so read natoms from TRX */
 992             if (!read_first_frame(oenv,&status,in_file,&fr,TRX_DONT_SKIP))
 993                 gmx_fatal(FARGS,"Could not read a frame from %s",in_file);
 994             natoms = fr.natoms;
 995             close_trj(status);
 996             sfree(fr.x);
 997             snew(index,natoms);
 998             for(i=0;i<natoms;i++)
 999                 index[i]=i;
1000             nout=natoms;
1001             if (bCenter) {
1002                 ncent = nout;
1003                 cindex = index;
1004             }
1005         }
1006
1007         if (bReset) {
1008             snew(w_rls,atoms->nr);
1009             for(i=0; (i<ifit); i++)
1010                 w_rls[ind_fit[i]]=atoms->atom[ind_fit[i]].m;
1011
1012             /* Restore reference structure and set to origin,
1013          store original location (to put structure back) */
1014             if (bRmPBC)
1015               gmx_rmpbc(gpbc,top.atoms.nr,top_box,xp);
1016             copy_rvec(xp[index[0]],x_shift);
1017             reset_x_ndim(nfitdim,ifit,ind_fit,atoms->nr,NULL,xp,w_rls);
1018             rvec_dec(x_shift,xp[index[0]]);
1019         } else
1020             clear_rvec(x_shift);
1021
1022         if (bDropUnder || bDropOver) {
1023             /* Read the xvg file with the drop values */
1024             fprintf(stderr,"\nReading drop file ...");
1025             ndrop = read_xvg(opt2fn("-drop",NFILE,fnm),&dropval,&ncol);
1026             fprintf(stderr," %d time points\n",ndrop);
1027             if (ndrop == 0 || ncol < 2)
1028                 gmx_fatal(FARGS,"Found no data points in %s",
1029                           opt2fn("-drop",NFILE,fnm));
1030             drop0 = 0;
1031             drop1 = 0;
1032         }
1033
1034         /* Make atoms struct for output in GRO or PDB files */
1035         if ((ftp == efGRO) || ((ftp == efG96) && bTPS) || (ftp == efPDB)) {
1036             /* get memory for stuff to go in pdb file */
1037             init_t_atoms(&useatoms,atoms->nr,FALSE);
1038             sfree(useatoms.resinfo);
1039             useatoms.resinfo = atoms->resinfo;
1040             for(i=0;(i<nout);i++) {
1041                 useatoms.atomname[i]=atoms->atomname[index[i]];
1042                 useatoms.atom[i]=atoms->atom[index[i]];
1043                 useatoms.nres=max(useatoms.nres,useatoms.atom[i].resind+1);
1044             }
1045             useatoms.nr=nout;
1046         }
1047         /* select what to read */
1048         if (ftp==efTRR || ftp==efTRJ)
1049             flags = TRX_READ_X;
1050         else
1051             flags = TRX_NEED_X;
1052         if (bVels)
1053             flags = flags | TRX_READ_V;
1054         if (bForce)
1055             flags = flags | TRX_READ_F;
1056
1057         /* open trx file for reading */
1058         bHaveFirstFrame = read_first_frame(oenv,&status,in_file,&fr,flags);
1059         if (fr.bPrec)
1060             fprintf(stderr,"\nPrecision of %s is %g (nm)\n",in_file,1/fr.prec);
1061         if (bNeedPrec) {
1062             if (bSetPrec || !fr.bPrec)
1063                 fprintf(stderr,"\nSetting output precision to %g (nm)\n",1/prec);
1064             else
1065                 fprintf(stderr,"Using output precision of %g (nm)\n",1/prec);
1066         }
1067
1068         if (bHaveFirstFrame) {
1069             set_trxframe_ePBC(&fr,ePBC);
1070
1071             natoms = fr.natoms;
1072
1073             if (bSetTime)
1074                 tshift=tzero-fr.time;
1075             else
1076                 tzero=fr.time;
1077
1078             /* open output for writing */
1079             if ((bAppend) && (gmx_fexist(out_file))) {
1080                 strcpy(filemode,"a");
1081                 fprintf(stderr,"APPENDING to existing file %s\n",out_file);
1082             } else
1083                 strcpy(filemode,"w");
1084             switch (ftp) {
1085             case efXTC:
1086             case efG87:
1087             case efTRR:
1088             case efTRJ:
1089                 out=NULL;
1090                 if (!bSplit && !bSubTraj)
1091                     trxout = open_trx(out_file,filemode);
1092                 break;
1093             case efGRO:
1094             case efG96:
1095             case efPDB:
1096                 if (( !bSeparate && !bSplit ) && !bSubTraj)
1097                     out=ffopen(out_file,filemode);
1098                 break;
1099             }
1100
1101             bCopy = FALSE;
1102             if (bIndex)
1103                 /* check if index is meaningful */
1104                 for(i=0; i<nout; i++) {
1105                     if (index[i] >= natoms)
1106                         gmx_fatal(FARGS,"Index[%d] %d is larger than the number of atoms in the trajectory file (%d)",i,index[i]+1,natoms);
1107                     bCopy = bCopy || (i != index[i]);
1108                 }
1109             if (bCopy) {
1110                 snew(xmem,nout);
1111                 if (bVels) {
1112                     snew(vmem,nout);
1113                 }
1114                 if (bForce) {
1115                     snew(fmem,nout);
1116                 }
1117             }
1118
1119             if (ftp == efG87)
1120                 fprintf(gmx_fio_getfp(trx_get_fileio(trxout)),
1121                         "Generated by %s. #atoms=%d, a BOX is"
1122                         " stored in this file.\n",Program(),nout);
1123
1124             /* Start the big loop over frames */
1125             file_nr  =  0;
1126             frame    =  0;
1127             outframe =  0;
1128             model_nr = -1;
1129             bDTset   = FALSE;
1130
1131             /* Main loop over frames */
1132             do {
1133                 if (!fr.bStep) {
1134                     /* set the step */
1135                     fr.step = newstep;
1136                     newstep++;
1137                 }
1138                 if (bSubTraj) {
1139                     /*if (frame >= clust->clust->nra)
1140             gmx_fatal(FARGS,"There are more frames in the trajectory than in the cluster index file\n");*/
1141                     if (frame >= clust->maxframe)
1142                         my_clust = -1;
1143                     else
1144                         my_clust = clust->inv_clust[frame];
1145                     if ((my_clust < 0) || (my_clust >= clust->clust->nr) ||
1146                         (my_clust == NO_ATID))
1147                         my_clust = -1;
1148                 }
1149
1150                 if (bSetBox) {
1151                     /* generate new box */
1152                     clear_mat(fr.box);
1153                     for (m=0; m<DIM; m++)
1154                         fr.box[m][m] = newbox[m];
1155                 }
1156
1157                 if (bTrans) {
1158                     for(i=0; i<natoms; i++)
1159                         rvec_inc(fr.x[i],trans);
1160                 }
1161
1162                 if (bTDump) {
1163                     /* determine timestep */
1164                     if (t0 == -1) {
1165                         t0 = fr.time;
1166                     } else {
1167                         if (!bDTset) {
1168                             dt = fr.time-t0;
1169                             bDTset = TRUE;
1170                         }
1171                     }
1172                     /* This is not very elegant, as one can not dump a frame after
1173                      * a timestep with is more than twice as small as the first one. */
1174                     bDumpFrame = (fr.time > tdump-0.5*dt) && (fr.time <= tdump+0.5*dt);
1175                 } else
1176                     bDumpFrame = FALSE;
1177
1178                 /* determine if an atom jumped across the box and reset it if so */
1179                 if (bNoJump && (bTPS || frame!=0)) {
1180                     for(d=0; d<DIM; d++)
1181                         hbox[d] = 0.5*fr.box[d][d];
1182                     for(i=0; i<natoms; i++) {
1183                         if (bReset)
1184                             rvec_dec(fr.x[i],x_shift);
1185                         for(m=DIM-1; m>=0; m--)
1186                             if (hbox[m] > 0) {
1187                                 while (fr.x[i][m]-xp[i][m] <= -hbox[m])
1188                                     for(d=0; d<=m; d++)
1189                                         fr.x[i][d] += fr.box[m][d];
1190                                 while (fr.x[i][m]-xp[i][m] > hbox[m])
1191                                     for(d=0; d<=m; d++)
1192                                         fr.x[i][d] -= fr.box[m][d];
1193                             }
1194                     }
1195                 }
1196                 else if (bCluster) {
1197                     rvec com;
1198
1199                     calc_pbc_cluster(ecenter,ifit,&top,ePBC,fr.x,ind_fit,com,fr.box);
1200                 }
1201
1202                 if (bPFit) {
1203                     /* Now modify the coords according to the flags,
1204              for normal fit, this is only done for output frames */
1205                     if (bRmPBC)
1206                       gmx_rmpbc_trxfr(gpbc,&fr);
1207
1208                     reset_x_ndim(nfitdim,ifit,ind_fit,natoms,NULL,fr.x,w_rls);
1209                     do_fit(natoms,w_rls,xp,fr.x);
1210                 }
1211
1212                 /* store this set of coordinates for future use */
1213                 if (bPFit || bNoJump) {
1214                     if (xp == NULL)
1215                         snew(xp,natoms);
1216                     for(i=0; (i<natoms); i++) {
1217                         copy_rvec(fr.x[i],xp[i]);
1218                         rvec_inc(fr.x[i],x_shift);
1219                     }
1220                 }
1221
1222                 if ( frindex )
1223                     /* see if we have a frame from the frame index group */
1224                     for(i=0; i<nrfri && !bDumpFrame; i++)
1225                         bDumpFrame = frame == frindex[i];
1226                 if (debug && bDumpFrame)
1227                     fprintf(debug,"dumping %d\n",frame);
1228
1229                 bWriteFrame =
1230                     ( ( !bTDump && !frindex && frame % skip_nr == 0 ) || bDumpFrame );
1231
1232                 if (bWriteFrame && (bDropUnder || bDropOver)) {
1233                     while (dropval[0][drop1]<fr.time && drop1+1<ndrop) {
1234                         drop0 = drop1;
1235                         drop1++;
1236                     }
1237                     if (fabs(dropval[0][drop0] - fr.time)
1238                         < fabs(dropval[0][drop1] - fr.time)) {
1239                         dropuse = drop0;
1240                     } else {
1241                         dropuse = drop1;
1242                     }
1243                     if ((bDropUnder && dropval[1][dropuse] < dropunder) ||
1244                         (bDropOver && dropval[1][dropuse] > dropover))
1245                         bWriteFrame = FALSE;
1246                 }
1247
1248                 if (bWriteFrame) {
1249
1250                     /* calc new time */
1251                     if (bTimeStep)
1252                         fr.time = tzero+frame*timestep;
1253                     else
1254                         if (bSetTime)
1255                             fr.time += tshift;
1256
1257                     if (bTDump)
1258                         fprintf(stderr,"\nDumping frame at t= %g %s\n",
1259                                 output_env_conv_time(oenv,fr.time),output_env_get_time_unit(oenv));
1260
1261                     /* check for writing at each delta_t */
1262                     bDoIt=(delta_t == 0);
1263                     if (!bDoIt)
1264                     {
1265                         if (!bRound)
1266                             bDoIt=bRmod(fr.time,tzero, delta_t);
1267                         else
1268                             /* round() is not C89 compatible, so we do this:  */
1269                             bDoIt=bRmod(floor(fr.time+0.5),floor(tzero+0.5),
1270                                         floor(delta_t+0.5));
1271                     }
1272
1273                     if (bDoIt || bTDump) {
1274                         /* print sometimes */
1275                         if ( ((outframe % SKIP) == 0) || (outframe < SKIP) )
1276                             fprintf(stderr," ->  frame %6d time %8.3f      \r",
1277                                     outframe,output_env_conv_time(oenv,fr.time));
1278
1279                         if (!bPFit) {
1280                             /* Now modify the coords according to the flags,
1281                                for PFit we did this already! */
1282
1283                             if (bRmPBC)
1284                               gmx_rmpbc_trxfr(gpbc,&fr);
1285
1286                             if (bReset) {
1287                                 reset_x_ndim(nfitdim,ifit,ind_fit,natoms,NULL,fr.x,w_rls);
1288                                 if (bFit)
1289                                     do_fit_ndim(nfitdim,natoms,w_rls,xp,fr.x);
1290                                 if (!bCenter)
1291                                     for(i=0; i<natoms; i++)
1292                                         rvec_inc(fr.x[i],x_shift);
1293                             }
1294
1295                             if (bCenter)
1296                                 center_x(ecenter,fr.x,fr.box,natoms,ncent,cindex);
1297                         }
1298
1299                         if (bPBCcomAtom) {
1300                             switch (unitcell_enum) {
1301                             case euRect:
1302                                 put_atoms_in_box(fr.box,natoms,fr.x);
1303                                 break;
1304                             case euTric:
1305                                 put_atoms_in_triclinic_unitcell(ecenter,fr.box,natoms,fr.x);
1306                                 break;
1307                             case euCompact:
1308                                 warn = put_atoms_in_compact_unitcell(ePBC,ecenter,fr.box,
1309                                                                      natoms,fr.x);
1310                                 if (warn && !bWarnCompact) {
1311                                     fprintf(stderr,"\n%s\n",warn);
1312                                     bWarnCompact = TRUE;
1313                                 }
1314                                 break;
1315                             }
1316                         }
1317                         if (bPBCcomRes) {
1318                             put_residue_com_in_box(unitcell_enum,ecenter,
1319                                                    natoms,atoms->atom,ePBC,fr.box,fr.x);
1320                         }
1321                         if (bPBCcomMol) {
1322                             put_molecule_com_in_box(unitcell_enum,ecenter,
1323                                                     &top.mols,
1324                                                     natoms,atoms->atom,ePBC,fr.box,fr.x);
1325                         }
1326                         /* Copy the input trxframe struct to the output trxframe struct */
1327                         frout = fr;
1328                         frout.natoms = nout;
1329                         if (bNeedPrec && (bSetPrec || !fr.bPrec)) {
1330                             frout.bPrec = TRUE;
1331                             frout.prec  = prec;
1332                         }
1333                         if (bCopy) {
1334                             frout.x = xmem;
1335                             if (bVels) {
1336                                 frout.v = vmem;
1337                             }
1338                             if (bForce) {
1339                                 frout.f = fmem;
1340                             }
1341                             for(i=0; i<nout; i++) {
1342                                 copy_rvec(fr.x[index[i]],frout.x[i]);
1343                                 if (bVels && fr.bV) {
1344                                     copy_rvec(fr.v[index[i]],frout.v[i]);
1345                                 }
1346                                 if (bForce && fr.bF) {
1347                                     copy_rvec(fr.f[index[i]],frout.f[i]);
1348                                 }
1349                             }
1350                         }
1351
1352                         if (opt2parg_bSet("-shift",NPA,pa))
1353                             for(i=0; i<nout; i++)
1354                                 for (d=0; d<DIM; d++)
1355                                     frout.x[i][d] += outframe*shift[d];
1356
1357                         if (!bRound)
1358                             bSplitHere = bSplit && bRmod(fr.time,tzero, split_t);
1359                         else
1360                         {
1361                             /* round() is not C89 compatible, so we do this: */
1362                             bSplitHere = bSplit && bRmod(floor(fr.time+0.5),
1363                                                          floor(tzero+0.5),
1364                                                          floor(split_t+0.5));
1365                         }
1366                         if (bSeparate || bSplitHere)
1367                             mk_filenm(outf_base,ftp2ext(ftp),nzero,file_nr,out_file2);
1368
1369                         switch(ftp) {
1370                         case efTRJ:
1371                         case efTRR:
1372                         case efG87:
1373                         case efXTC:
1374                             if ( bSplitHere ) {
1375                                 if ( trxout )
1376                                     close_trx(trxout);
1377                                 trxout = open_trx(out_file2,filemode);
1378                             }
1379                             if (bSubTraj) {
1380                                 if (my_clust != -1) {
1381                                     char buf[STRLEN];
1382                                     if (clust_status_id[my_clust] == -1) {
1383                                         sprintf(buf,"%s.%s",clust->grpname[my_clust],ftp2ext(ftp));
1384                                         clust_status[my_clust] = open_trx(buf,"w");
1385                                         clust_status_id[my_clust] = 1;
1386                                         ntrxopen++;
1387                                     }
1388                                     else if (clust_status_id[my_clust] == -2)
1389                                         gmx_fatal(FARGS,"File %s.xtc should still be open (%d open xtc files)\n""in order to write frame %d. my_clust = %d",
1390                                                   clust->grpname[my_clust],ntrxopen,frame,
1391                                                   my_clust);
1392                                     write_trxframe(clust_status[my_clust],&frout,gc);
1393                                     nfwritten[my_clust]++;
1394                                     if (nfwritten[my_clust] ==
1395                                         (clust->clust->index[my_clust+1]-
1396                                             clust->clust->index[my_clust])) {
1397                                         close_trx(clust_status[my_clust]);
1398                                         clust_status_id[my_clust] = -2;
1399                                         ntrxopen--;
1400                                         if (ntrxopen < 0)
1401                                             gmx_fatal(FARGS,"Less than zero open xtc files!");
1402                                     }
1403                                 }
1404                             }
1405                             else
1406                                 write_trxframe(trxout,&frout,gc);
1407                             break;
1408                         case efGRO:
1409                         case efG96:
1410                         case efPDB:
1411                             sprintf(title,"Generated by trjconv : %s t= %9.5f",
1412                                     top_title,fr.time);
1413                             if (bSeparate || bSplitHere)
1414                                 out=ffopen(out_file2,"w");
1415                             switch(ftp) {
1416                             case efGRO:
1417                                 write_hconf_p(out,title,&useatoms,prec2ndec(frout.prec),
1418                                               frout.x,fr.bV?frout.v:NULL,frout.box);
1419                                 break;
1420                             case efPDB:
1421                                 fprintf(out,"REMARK    GENERATED BY TRJCONV\n");
1422                                 sprintf(title,"%s t= %9.5f",top_title,frout.time);
1423                                 /* if reading from pdb, we want to keep the original
1424                    model numbering else we write the output frame
1425                    number plus one, because model 0 is not allowed in pdb */
1426                                 if (ftpin==efPDB && fr.bStep && fr.step > model_nr)
1427                                     model_nr = fr.step;
1428                                 else
1429                                     model_nr++;
1430                                 write_pdbfile(out,title,&useatoms,frout.x,
1431                                               frout.ePBC,frout.box,' ',model_nr,gc,TRUE);
1432                                 break;
1433                             case efG96:
1434                                 frout.title = title;
1435                                 if (bSeparate || bTDump) {
1436                                     frout.bTitle = TRUE;
1437                                     if (bTPS)
1438                                         frout.bAtoms = TRUE;
1439                                     frout.atoms  = &useatoms;
1440                                     frout.bStep  = FALSE;
1441                                     frout.bTime  = FALSE;
1442                                 } else {
1443                                     frout.bTitle = (outframe == 0);
1444                                     frout.bAtoms = FALSE;
1445                                     frout.bStep = TRUE;
1446                                     frout.bTime = TRUE;
1447                                 }
1448                                 write_g96_conf(out,&frout,-1,NULL);
1449                             }
1450                             if (bSeparate) {
1451                                 ffclose(out);
1452                                 out = NULL;
1453                             }
1454                             break;
1455                             default:
1456                                 gmx_fatal(FARGS,"DHE, ftp=%d\n",ftp);
1457                         }
1458                         if (bSeparate || bSplitHere)
1459                             file_nr++;
1460
1461                         /* execute command */
1462                         if (bExec) {
1463                             char c[255];
1464                             sprintf(c,"%s  %d",exec_command,file_nr-1);
1465                             /*fprintf(stderr,"Executing '%s'\n",c);*/
1466 #ifdef GMX_NO_SYSTEM
1467                             printf("Warning-- No calls to system(3) supported on this platform.");
1468                             printf("Warning-- Skipping execution of 'system(\"%s\")'.", c);
1469 #else
1470                             if(0 != system(c))
1471                             {
1472                                 gmx_fatal(FARGS,"Error executing command: %s",c);
1473                             }
1474 #endif
1475                         }
1476                         outframe++;
1477                     }
1478                 }
1479                 frame++;
1480                 bHaveNextFrame=read_next_frame(oenv,status,&fr);
1481             } while (!(bTDump && bDumpFrame) && bHaveNextFrame);
1482         }
1483
1484         if (!bHaveFirstFrame || (bTDump && !bDumpFrame))
1485             fprintf(stderr,"\nWARNING no output, "
1486                     "last frame read at t=%g\n",fr.time);
1487         fprintf(stderr,"\n");
1488
1489         close_trj(status);
1490         if (bRmPBC)
1491           gmx_rmpbc_done(gpbc);
1492
1493         if (trxout)
1494             close_trx(trxout);
1495         else if (out != NULL)
1496             ffclose(out);
1497         if (bSubTraj) {
1498             for(i=0; (i<clust->clust->nr); i++)
1499                 if (clust_status[i] )
1500                     close_trx(clust_status[i]);
1501         }
1502     }
1503
1504     do_view(oenv,out_file,NULL);
1505
1506     thanx(stderr);
1507
1508     return 0;
1509 }