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