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