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(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(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     CopyRight(stderr,argv[0]);
 818     parse_common_args(&argc,argv,
 819                       PCA_CAN_BEGIN | PCA_CAN_END | PCA_CAN_VIEW |
 820                       PCA_TIME_UNIT | PCA_BE_NICE,
 821                       NFILE,fnm,NPA,pa,asize(desc),desc,
 822                       0,NULL,&oenv);
 823
 824     top_file = ftp2fn(efTPS,NFILE,fnm);
 825     init_top(&top);
 826
 827     /* Check command line */
 828     in_file=opt2fn("-f",NFILE,fnm);
 829
 830     if (ttrunc != -1) {
 831 #ifndef GMX_NATIVE_WINDOWS
 832         do_trunc(in_file,ttrunc);
 833 #endif
 834     }
 835     else {
 836         /* mark active cmdline options */
 837         bSetBox   = opt2parg_bSet("-box", NPA, pa);
 838         bSetTime  = opt2parg_bSet("-t0", NPA, pa);
 839         bSetPrec  = opt2parg_bSet("-ndec", NPA, pa);
 840         bSetUR    = opt2parg_bSet("-ur", NPA, pa);
 841         bExec     = opt2parg_bSet("-exec", NPA, pa);
 842         bTimeStep = opt2parg_bSet("-timestep", NPA, pa);
 843         bTDump    = opt2parg_bSet("-dump", NPA, pa);
 844         bDropUnder= opt2parg_bSet("-dropunder", NPA, pa);
 845         bDropOver = opt2parg_bSet("-dropover", NPA, pa);
 846         bTrans    = opt2parg_bSet("-trans",NPA,pa);
 847         bSplit    = (split_t != 0);
 848
 849         /* parse enum options */
 850         fit_enum   = nenum(fit);
 851         bFit       = (fit_enum==efFit || fit_enum==efFitXY);
 852         bFitXY     = fit_enum==efFitXY;
 853         bReset     = (fit_enum==efReset || fit_enum==efResetXY);
 854         bPFit      = fit_enum==efPFit;
 855         pbc_enum   = nenum(pbc_opt);
 856         bPBCWhole  = pbc_enum==epWhole;
 857         bPBCcomRes = pbc_enum==epComRes;
 858         bPBCcomMol = pbc_enum==epComMol;
 859         bPBCcomAtom= pbc_enum==epComAtom;
 860         bNoJump    = pbc_enum==epNojump;
 861         bCluster   = pbc_enum==epCluster;
 862         bPBC       = pbc_enum!=epNone;
 863         unitcell_enum = nenum(unitcell_opt);
 864         ecenter    = nenum(center_opt) - ecTric;
 865
 866         /* set and check option dependencies */
 867         if (bPFit) bFit = TRUE; /* for pfit, fit *must* be set */
 868         if (bFit) bReset = TRUE; /* for fit, reset *must* be set */
 869         nfitdim = 0;
 870         if (bFit || bReset)
 871             nfitdim = (fit_enum==efFitXY || fit_enum==efResetXY) ? 2 : 3;
 872         bRmPBC = bFit || bPBCWhole || bPBCcomRes || bPBCcomMol;
 873
 874         if (bSetUR) {
 875             if (!(bPBCcomRes || bPBCcomMol ||  bPBCcomAtom)) {
 876                 fprintf(stderr,
 877                         "WARNING: Option for unitcell representation (-ur %s)\n"
 878                         "         only has effect in combination with -pbc %s, %s or %s.\n"
 879                         "         Ingoring unitcell representation.\n\n",
 880                         unitcell_opt[0],pbc_opt[2],pbc_opt[3],pbc_opt[4]);
 881                 bSetUR = FALSE;
 882             }
 883         }
 884         if (bFit && bPBC) {
 885             gmx_fatal(FARGS,"PBC condition treatment does not work together with rotational fit.\n"
 886                       "Please do the PBC condition treatment first and then run trjconv in a second step\n"
 887                       "for the rotational fit.\n"
 888                       "First doing the rotational fit and then doing the PBC treatment gives incorrect\n"
 889                       "results!");
 890         }
 891
 892         /* ndec is in nr of decimal places, prec is a multiplication factor: */
 893         prec = 1;
 894         for (i=0; i<ndec; i++)
 895             prec *= 10;
 896
 897         bIndex=ftp2bSet(efNDX,NFILE,fnm);
 898
 899
 900         /* Determine output type */
 901         out_file=opt2fn("-o",NFILE,fnm);
 902         ftp=fn2ftp(out_file);
 903         fprintf(stderr,"Will write %s: %s\n",ftp2ext(ftp),ftp2desc(ftp));
 904         bNeedPrec = (ftp==efXTC || ftp==efGRO);
 905         if (bVels) {
 906             /* check if velocities are possible in input and output files */
 907             ftpin=fn2ftp(in_file);
 908             bVels= (ftp==efTRR || ftp==efTRJ || ftp==efGRO || ftp==efG96)
 909                 && (ftpin==efTRR || ftpin==efTRJ || ftpin==efGRO || ftpin==efG96 ||
 910                     ftpin==efCPT);
 911         }
 912         if (bSeparate || bSplit) {
 913             outf_ext = strrchr(out_file,'.');
 914             if (outf_ext == NULL)
 915                 gmx_fatal(FARGS,"Output file name '%s' does not contain a '.'",out_file);
 916             outf_base = strdup(out_file);
 917             outf_base[outf_ext - out_file] = '\0';
 918         }
 919
 920         bSubTraj = opt2bSet("-sub",NFILE,fnm);
 921         if (bSubTraj) {
 922             if ((ftp != efXTC) && (ftp != efTRN))
 923                 gmx_fatal(FARGS,"Can only use the sub option with output file types "
 924                           "xtc and trr");
 925             clust = cluster_index(NULL,opt2fn("-sub",NFILE,fnm));
 926
 927             /* Check for number of files disabled, as FOPEN_MAX is not the correct
 928              * number to check for. In my linux box it is only 16.
 929              */
 930             if (0 && (clust->clust->nr > FOPEN_MAX-4))
 931                 gmx_fatal(FARGS,"Can not open enough (%d) files to write all the"
 932                           " trajectories.\ntry splitting the index file in %d parts.\n"
 933                           "FOPEN_MAX = %d",
 934                           clust->clust->nr,1+clust->clust->nr/FOPEN_MAX,FOPEN_MAX);
 935             gmx_warning("The -sub option could require as many open output files as there are\n"
 936                         "index groups in the file (%d). If you get I/O errors opening new files,\n"
 937                         "try reducing the number of index groups in the file, and perhaps\n"
 938                         "using trjconv -sub several times on different chunks of your index file.\n",
 939                         clust->clust->nr);
 940
 941             snew(clust_status,clust->clust->nr);
 942             snew(clust_status_id,clust->clust->nr);
 943             snew(nfwritten,clust->clust->nr);
 944             for(i=0; (i<clust->clust->nr); i++)
 945             {
 946                 clust_status[i] = NULL;
 947                 clust_status_id[i] = -1;
 948             }
 949             bSeparate = bSplit = FALSE;
 950         }
 951         /* skipping */
 952         if (skip_nr <= 0) {
 953         }
 954
 955         /* Determine whether to read a topology */
 956         bTPS = (ftp2bSet(efTPS,NFILE,fnm) ||
 957             bRmPBC || bReset || bPBCcomMol || bCluster ||
 958             (ftp == efGRO) || (ftp == efPDB) || bCONECT);
 959
 960         /* Determine if when can read index groups */
 961         bIndex = (bIndex || bTPS);
 962
 963         if (bTPS) {
 964             read_tps_conf(top_file,top_title,&top,&ePBC,&xp,NULL,top_box,
 965                           bReset || bPBCcomRes);
 966             atoms=&top.atoms;
 967
 968             if (0 == top.mols.nr && (bCluster || bPBCcomMol))
 969             {
 970                 gmx_fatal(FARGS,"Option -pbc %s requires a .tpr file for the -s option",pbc_opt[pbc_enum]);
 971             }
 972
 973             /* top_title is only used for gro and pdb,
 974              * the header in such a file is top_title t= ...
 975              * to prevent a double t=, remove it from top_title
 976              */
 977             if ((charpt=strstr(top_title," t= ")))
 978                 charpt[0]='\0';
 979
 980             if (bCONECT)
 981                 gc = gmx_conect_generate(&top);
 982             if (bRmPBC)
 983               gpbc = gmx_rmpbc_init(&top.idef,ePBC,top.atoms.nr,top_box);
 984         }
 985
 986         /* get frame number index */
 987         frindex=NULL;
 988         if (opt2bSet("-fr",NFILE,fnm)) {
 989             printf("Select groups of frame number indices:\n");
 990             rd_index(opt2fn("-fr",NFILE,fnm),1,&nrfri,(atom_id **)&frindex,&frname);
 991             if (debug)
 992                 for(i=0; i<nrfri; i++)
 993                     fprintf(debug,"frindex[%4d]=%4d\n",i,frindex[i]);
 994         }
 995
 996         /* get index groups etc. */
 997         if (bReset) {
 998             printf("Select group for %s fit\n",
 999                    bFit?"least squares":"translational");
1000             get_index(atoms,ftp2fn_null(efNDX,NFILE,fnm),
1001                       1,&ifit,&ind_fit,&gn_fit);
1002
1003             if (bFit) {
1004                 if (ifit < 2)
1005                     gmx_fatal(FARGS,"Need at least 2 atoms to fit!\n");
1006                 else if (ifit == 3)
1007                     fprintf(stderr,"WARNING: fitting with only 2 atoms is not unique\n");
1008             }
1009         }
1010         else if (bCluster) {
1011             printf("Select group for clustering\n");
1012             get_index(atoms,ftp2fn_null(efNDX,NFILE,fnm),
1013                       1,&ifit,&ind_fit,&gn_fit);
1014         }
1015
1016         if (bIndex) {
1017             if (bCenter) {
1018                 printf("Select group for centering\n");
1019                 get_index(atoms,ftp2fn_null(efNDX,NFILE,fnm),
1020                           1,&ncent,&cindex,&grpnm);
1021             }
1022             printf("Select group for output\n");
1023             get_index(atoms,ftp2fn_null(efNDX,NFILE,fnm),
1024                       1,&nout,&index,&grpnm);
1025         } else {
1026             /* no index file, so read natoms from TRX */
1027             if (!read_first_frame(oenv,&status,in_file,&fr,TRX_DONT_SKIP))
1028                 gmx_fatal(FARGS,"Could not read a frame from %s",in_file);
1029             natoms = fr.natoms;
1030             close_trj(status);
1031             sfree(fr.x);
1032             snew(index,natoms);
1033             for(i=0;i<natoms;i++)
1034                 index[i]=i;
1035             nout=natoms;
1036             if (bCenter) {
1037                 ncent = nout;
1038                 cindex = index;
1039             }
1040         }
1041
1042         if (bReset) {
1043             snew(w_rls,atoms->nr);
1044             for(i=0; (i<ifit); i++)
1045                 w_rls[ind_fit[i]]=atoms->atom[ind_fit[i]].m;
1046
1047             /* Restore reference structure and set to origin,
1048          store original location (to put structure back) */
1049             if (bRmPBC)
1050               gmx_rmpbc(gpbc,top.atoms.nr,top_box,xp);
1051             copy_rvec(xp[index[0]],x_shift);
1052             reset_x_ndim(nfitdim,ifit,ind_fit,atoms->nr,NULL,xp,w_rls);
1053             rvec_dec(x_shift,xp[index[0]]);
1054         } else
1055             clear_rvec(x_shift);
1056
1057         if (bDropUnder || bDropOver) {
1058             /* Read the .xvg file with the drop values */
1059             fprintf(stderr,"\nReading drop file ...");
1060             ndrop = read_xvg(opt2fn("-drop",NFILE,fnm),&dropval,&ncol);
1061             fprintf(stderr," %d time points\n",ndrop);
1062             if (ndrop == 0 || ncol < 2)
1063                 gmx_fatal(FARGS,"Found no data points in %s",
1064                           opt2fn("-drop",NFILE,fnm));
1065             drop0 = 0;
1066             drop1 = 0;
1067         }
1068
1069         /* Make atoms struct for output in GRO or PDB files */
1070         if ((ftp == efGRO) || ((ftp == efG96) && bTPS) || (ftp == efPDB)) {
1071             /* get memory for stuff to go in .pdb file */
1072             init_t_atoms(&useatoms,atoms->nr,FALSE);
1073             sfree(useatoms.resinfo);
1074             useatoms.resinfo = atoms->resinfo;
1075             for(i=0;(i<nout);i++) {
1076                 useatoms.atomname[i]=atoms->atomname[index[i]];
1077                 useatoms.atom[i]=atoms->atom[index[i]];
1078                 useatoms.nres=max(useatoms.nres,useatoms.atom[i].resind+1);
1079             }
1080             useatoms.nr=nout;
1081         }
1082         /* select what to read */
1083         if (ftp==efTRR || ftp==efTRJ)
1084             flags = TRX_READ_X;
1085         else
1086             flags = TRX_NEED_X;
1087         if (bVels)
1088             flags = flags | TRX_READ_V;
1089         if (bForce)
1090             flags = flags | TRX_READ_F;
1091
1092         /* open trx file for reading */
1093         bHaveFirstFrame = read_first_frame(oenv,&status,in_file,&fr,flags);
1094         if (fr.bPrec)
1095             fprintf(stderr,"\nPrecision of %s is %g (nm)\n",in_file,1/fr.prec);
1096         if (bNeedPrec) {
1097             if (bSetPrec || !fr.bPrec)
1098                 fprintf(stderr,"\nSetting output precision to %g (nm)\n",1/prec);
1099             else
1100                 fprintf(stderr,"Using output precision of %g (nm)\n",1/prec);
1101         }
1102
1103         if (bHaveFirstFrame) {
1104             set_trxframe_ePBC(&fr,ePBC);
1105
1106             natoms = fr.natoms;
1107
1108             if (bSetTime)
1109                 tshift=tzero-fr.time;
1110             else
1111                 tzero=fr.time;
1112
1113             /* open output for writing */
1114             if ((bAppend) && (gmx_fexist(out_file))) {
1115                 strcpy(filemode,"a");
1116                 fprintf(stderr,"APPENDING to existing file %s\n",out_file);
1117             } else
1118                 strcpy(filemode,"w");
1119             switch (ftp) {
1120             case efXTC:
1121             case efG87:
1122             case efTRR:
1123             case efTRJ:
1124                 out=NULL;
1125                 if (!bSplit && !bSubTraj)
1126                     trxout = open_trx(out_file,filemode);
1127                 break;
1128             case efGRO:
1129             case efG96:
1130             case efPDB:
1131                 if (( !bSeparate && !bSplit ) && !bSubTraj)
1132                     out=ffopen(out_file,filemode);
1133                 break;
1134             }
1135
1136             bCopy = FALSE;
1137             if (bIndex)
1138                 /* check if index is meaningful */
1139                 for(i=0; i<nout; i++) {
1140                     if (index[i] >= natoms)
1141                         gmx_fatal(FARGS,
1142                                   "Index[%d] %d is larger than the number of atoms in the\n"
1143                                   "trajectory file (%d). There is a mismatch in the contents\n"
1144                                   "of your -f, -s and/or -n files.",i,index[i]+1,natoms);
1145                     bCopy = bCopy || (i != index[i]);
1146                 }
1147             if (bCopy) {
1148                 snew(xmem,nout);
1149                 if (bVels) {
1150                     snew(vmem,nout);
1151                 }
1152                 if (bForce) {
1153                     snew(fmem,nout);
1154                 }
1155             }
1156
1157             if (ftp == efG87)
1158                 fprintf(gmx_fio_getfp(trx_get_fileio(trxout)),
1159                         "Generated by %s. #atoms=%d, a BOX is"
1160                         " stored in this file.\n",Program(),nout);
1161
1162             /* Start the big loop over frames */
1163             file_nr  =  0;
1164             frame    =  0;
1165             outframe =  0;
1166             model_nr =  0;
1167             bDTset   = FALSE;
1168
1169             /* Main loop over frames */
1170             do {
1171                 if (!fr.bStep) {
1172                     /* set the step */
1173                     fr.step = newstep;
1174                     newstep++;
1175                 }
1176                 if (bSubTraj) {
1177                     /*if (frame >= clust->clust->nra)
1178             gmx_fatal(FARGS,"There are more frames in the trajectory than in the cluster index file\n");*/
1179                     if (frame > clust->maxframe)
1180                         my_clust = -1;
1181                     else
1182                         my_clust = clust->inv_clust[frame];
1183                     if ((my_clust < 0) || (my_clust >= clust->clust->nr) ||
1184                         (my_clust == NO_ATID))
1185                         my_clust = -1;
1186                 }
1187
1188                 if (bSetBox) {
1189                     /* generate new box */
1190                     clear_mat(fr.box);
1191                     for (m=0; m<DIM; m++)
1192                         fr.box[m][m] = newbox[m];
1193                 }
1194
1195                 if (bTrans) {
1196                     for(i=0; i<natoms; i++)
1197                         rvec_inc(fr.x[i],trans);
1198                 }
1199
1200                 if (bTDump) {
1201                     /* determine timestep */
1202                     if (t0 == -1) {
1203                         t0 = fr.time;
1204                     } else {
1205                         if (!bDTset) {
1206                             dt = fr.time-t0;
1207                             bDTset = TRUE;
1208                         }
1209                     }
1210                     /* This is not very elegant, as one can not dump a frame after
1211                      * a timestep with is more than twice as small as the first one. */
1212                     bDumpFrame = (fr.time > tdump-0.5*dt) && (fr.time <= tdump+0.5*dt);
1213                 } else
1214                     bDumpFrame = FALSE;
1215
1216                 /* determine if an atom jumped across the box and reset it if so */
1217                 if (bNoJump && (bTPS || frame!=0)) {
1218                     for(d=0; d<DIM; d++)
1219                         hbox[d] = 0.5*fr.box[d][d];
1220                     for(i=0; i<natoms; i++) {
1221                         if (bReset)
1222                             rvec_dec(fr.x[i],x_shift);
1223                         for(m=DIM-1; m>=0; m--)
1224                             if (hbox[m] > 0) {
1225                                 while (fr.x[i][m]-xp[i][m] <= -hbox[m])
1226                                     for(d=0; d<=m; d++)
1227                                         fr.x[i][d] += fr.box[m][d];
1228                                 while (fr.x[i][m]-xp[i][m] > hbox[m])
1229                                     for(d=0; d<=m; d++)
1230                                         fr.x[i][d] -= fr.box[m][d];
1231                             }
1232                     }
1233                 }
1234                 else if (bCluster) {
1235                     rvec com;
1236
1237                     calc_pbc_cluster(ecenter,ifit,&top,ePBC,fr.x,ind_fit,com,fr.box,clustercenter);
1238                 }
1239
1240                 if (bPFit) {
1241                     /* Now modify the coords according to the flags,
1242              for normal fit, this is only done for output frames */
1243                     if (bRmPBC)
1244                       gmx_rmpbc_trxfr(gpbc,&fr);
1245
1246                     reset_x_ndim(nfitdim,ifit,ind_fit,natoms,NULL,fr.x,w_rls);
1247                     do_fit(natoms,w_rls,xp,fr.x);
1248                 }
1249
1250                 /* store this set of coordinates for future use */
1251                 if (bPFit || bNoJump) {
1252                     if (xp == NULL)
1253                         snew(xp,natoms);
1254                     for(i=0; (i<natoms); i++) {
1255                         copy_rvec(fr.x[i],xp[i]);
1256                         rvec_inc(fr.x[i],x_shift);
1257                     }
1258                 }
1259
1260                 if ( frindex )
1261                     /* see if we have a frame from the frame index group */
1262                     for(i=0; i<nrfri && !bDumpFrame; i++)
1263                         bDumpFrame = frame == frindex[i];
1264                 if (debug && bDumpFrame)
1265                     fprintf(debug,"dumping %d\n",frame);
1266
1267                 bWriteFrame =
1268                     ( ( !bTDump && !frindex && frame % skip_nr == 0 ) || bDumpFrame );
1269
1270                 if (bWriteFrame && (bDropUnder || bDropOver)) {
1271                     while (dropval[0][drop1]<fr.time && drop1+1<ndrop) {
1272                         drop0 = drop1;
1273                         drop1++;
1274                     }
1275                     if (fabs(dropval[0][drop0] - fr.time)
1276                         < fabs(dropval[0][drop1] - fr.time)) {
1277                         dropuse = drop0;
1278                     } else {
1279                         dropuse = drop1;
1280                     }
1281                     if ((bDropUnder && dropval[1][dropuse] < dropunder) ||
1282                         (bDropOver && dropval[1][dropuse] > dropover))
1283                         bWriteFrame = FALSE;
1284                 }
1285
1286                 if (bWriteFrame) {
1287
1288                     /* calc new time */
1289                     if (bTimeStep)
1290                         fr.time = tzero+frame*timestep;
1291                     else
1292                         if (bSetTime)
1293                             fr.time += tshift;
1294
1295                     if (bTDump)
1296                         fprintf(stderr,"\nDumping frame at t= %g %s\n",
1297                                 output_env_conv_time(oenv,fr.time),output_env_get_time_unit(oenv));
1298
1299                     /* check for writing at each delta_t */
1300                     bDoIt=(delta_t == 0);
1301                     if (!bDoIt)
1302                     {
1303                         if (!bRound)
1304                             bDoIt=bRmod(fr.time,tzero, delta_t);
1305                         else
1306                             /* round() is not C89 compatible, so we do this:  */
1307                             bDoIt=bRmod(floor(fr.time+0.5),floor(tzero+0.5),
1308                                         floor(delta_t+0.5));
1309                     }
1310
1311                     if (bDoIt || bTDump) {
1312                         /* print sometimes */
1313                         if ( ((outframe % SKIP) == 0) || (outframe < SKIP) )
1314                             fprintf(stderr," ->  frame %6d time %8.3f      \r",
1315                                     outframe,output_env_conv_time(oenv,fr.time));
1316
1317                         if (!bPFit) {
1318                             /* Now modify the coords according to the flags,
1319                                for PFit we did this already! */
1320
1321                             if (bRmPBC)
1322                               gmx_rmpbc_trxfr(gpbc,&fr);
1323
1324                             if (bReset) {
1325                                 reset_x_ndim(nfitdim,ifit,ind_fit,natoms,NULL,fr.x,w_rls);
1326                                 if (bFit)
1327                                     do_fit_ndim(nfitdim,natoms,w_rls,xp,fr.x);
1328                                 if (!bCenter)
1329                                     for(i=0; i<natoms; i++)
1330                                         rvec_inc(fr.x[i],x_shift);
1331                             }
1332
1333                             if (bCenter)
1334                                 center_x(ecenter,fr.x,fr.box,natoms,ncent,cindex);
1335                         }
1336
1337                         if (bPBCcomAtom) {
1338                             switch (unitcell_enum) {
1339                             case euRect:
1340                                 put_atoms_in_box(fr.box,natoms,fr.x);
1341                                 break;
1342                             case euTric:
1343                                 put_atoms_in_triclinic_unitcell(ecenter,fr.box,natoms,fr.x);
1344                                 break;
1345                             case euCompact:
1346                                 warn = put_atoms_in_compact_unitcell(ePBC,ecenter,fr.box,
1347                                                                      natoms,fr.x);
1348                                 if (warn && !bWarnCompact) {
1349                                     fprintf(stderr,"\n%s\n",warn);
1350                                     bWarnCompact = TRUE;
1351                                 }
1352                                 break;
1353                             }
1354                         }
1355                         if (bPBCcomRes) {
1356                             put_residue_com_in_box(unitcell_enum,ecenter,
1357                                                    natoms,atoms->atom,ePBC,fr.box,fr.x);
1358                         }
1359                         if (bPBCcomMol) {
1360                             put_molecule_com_in_box(unitcell_enum,ecenter,
1361                                                     &top.mols,
1362                                                     natoms,atoms->atom,ePBC,fr.box,fr.x);
1363                         }
1364                         /* Copy the input trxframe struct to the output trxframe struct */
1365                         frout = fr;
1366                         frout.bV = (frout.bV && bVels);
1367                         frout.bF = (frout.bF && bForce);
1368                         frout.natoms = nout;
1369                         if (bNeedPrec && (bSetPrec || !fr.bPrec)) {
1370                             frout.bPrec = TRUE;
1371                             frout.prec  = prec;
1372                         }
1373                         if (bCopy) {
1374                             frout.x = xmem;
1375                             if (frout.bV) {
1376                                 frout.v = vmem;
1377                             }
1378                             if (frout.bF) {
1379                                 frout.f = fmem;
1380                             }
1381                             for(i=0; i<nout; i++) {
1382                                 copy_rvec(fr.x[index[i]],frout.x[i]);
1383                                 if (frout.bV) {
1384                                     copy_rvec(fr.v[index[i]],frout.v[i]);
1385                                 }
1386                                 if (frout.bF) {
1387                                     copy_rvec(fr.f[index[i]],frout.f[i]);
1388                                 }
1389                             }
1390                         }
1391
1392                         if (opt2parg_bSet("-shift",NPA,pa))
1393                             for(i=0; i<nout; i++)
1394                                 for (d=0; d<DIM; d++)
1395                                     frout.x[i][d] += outframe*shift[d];
1396
1397                         if (!bRound)
1398                             bSplitHere = bSplit && bRmod(fr.time,tzero, split_t);
1399                         else
1400                         {
1401                             /* round() is not C89 compatible, so we do this: */
1402                             bSplitHere = bSplit && bRmod(floor(fr.time+0.5),
1403                                                          floor(tzero+0.5),
1404                                                          floor(split_t+0.5));
1405                         }
1406                         if (bSeparate || bSplitHere)
1407                             mk_filenm(outf_base,ftp2ext(ftp),nzero,file_nr,out_file2);
1408
1409                         switch(ftp) {
1410                         case efTRJ:
1411                         case efTRR:
1412                         case efG87:
1413                         case efXTC:
1414                             if ( bSplitHere ) {
1415                                 if ( trxout )
1416                                     close_trx(trxout);
1417                                 trxout = open_trx(out_file2,filemode);
1418                             }
1419                             if (bSubTraj) {
1420                                 if (my_clust != -1) {
1421                                     char buf[STRLEN];
1422                                     if (clust_status_id[my_clust] == -1) {
1423                                         sprintf(buf,"%s.%s",clust->grpname[my_clust],ftp2ext(ftp));
1424                                         clust_status[my_clust] = open_trx(buf,"w");
1425                                         clust_status_id[my_clust] = 1;
1426                                         ntrxopen++;
1427                                     }
1428                                     else if (clust_status_id[my_clust] == -2)
1429                                         gmx_fatal(FARGS,"File %s.xtc should still be open (%d open .xtc files)\n""in order to write frame %d. my_clust = %d",
1430                                                   clust->grpname[my_clust],ntrxopen,frame,
1431                                                   my_clust);
1432                                     write_trxframe(clust_status[my_clust],&frout,gc);
1433                                     nfwritten[my_clust]++;
1434                                     if (nfwritten[my_clust] ==
1435                                         (clust->clust->index[my_clust+1]-
1436                                             clust->clust->index[my_clust])) {
1437                                         close_trx(clust_status[my_clust]);
1438                                         clust_status[my_clust] = NULL;
1439                                         clust_status_id[my_clust] = -2;
1440                                         ntrxopen--;
1441                                         if (ntrxopen < 0)
1442                                             gmx_fatal(FARGS,"Less than zero open .xtc files!");
1443                                     }
1444                                 }
1445                             }
1446                             else
1447                                 write_trxframe(trxout,&frout,gc);
1448                             break;
1449                         case efGRO:
1450                         case efG96:
1451                         case efPDB:
1452                             sprintf(title,"Generated by trjconv : %s t= %9.5f",
1453                                     top_title,fr.time);
1454                             if (bSeparate || bSplitHere)
1455                                 out=ffopen(out_file2,"w");
1456                             switch(ftp) {
1457                             case efGRO:
1458                                 write_hconf_p(out,title,&useatoms,prec2ndec(frout.prec),
1459                                               frout.x,frout.bV?frout.v:NULL,frout.box);
1460                                 break;
1461                             case efPDB:
1462                                 fprintf(out,"REMARK    GENERATED BY TRJCONV\n");
1463                                 sprintf(title,"%s t= %9.5f",top_title,frout.time);
1464                                 /* if reading from pdb, we want to keep the original
1465                    model numbering else we write the output frame
1466                    number plus one, because model 0 is not allowed in pdb */
1467                                 if (ftpin==efPDB && fr.bStep && fr.step > model_nr)
1468                                     model_nr = fr.step;
1469                                 else
1470                                     model_nr++;
1471                                 write_pdbfile(out,title,&useatoms,frout.x,
1472                                               frout.ePBC,frout.box,' ',model_nr,gc,TRUE);
1473                                 break;
1474                             case efG96:
1475                                 frout.title = title;
1476                                 if (bSeparate || bTDump) {
1477                                     frout.bTitle = TRUE;
1478                                     if (bTPS)
1479                                         frout.bAtoms = TRUE;
1480                                     frout.atoms  = &useatoms;
1481                                     frout.bStep  = FALSE;
1482                                     frout.bTime  = FALSE;
1483                                 } else {
1484                                     frout.bTitle = (outframe == 0);
1485                                     frout.bAtoms = FALSE;
1486                                     frout.bStep = TRUE;
1487                                     frout.bTime = TRUE;
1488                                 }
1489                                 write_g96_conf(out,&frout,-1,NULL);
1490                             }
1491                             if (bSeparate) {
1492                                 ffclose(out);
1493                                 out = NULL;
1494                             }
1495                             break;
1496                             default:
1497                                 gmx_fatal(FARGS,"DHE, ftp=%d\n",ftp);
1498                         }
1499                         if (bSeparate || bSplitHere)
1500                             file_nr++;
1501
1502                         /* execute command */
1503                         if (bExec) {
1504                             char c[255];
1505                             sprintf(c,"%s  %d",exec_command,file_nr-1);
1506                             /*fprintf(stderr,"Executing '%s'\n",c);*/
1507 #ifdef GMX_NO_SYSTEM
1508                             printf("Warning-- No calls to system(3) supported on this platform.");
1509                             printf("Warning-- Skipping execution of 'system(\"%s\")'.", c);
1510 #else
1511                             if(0 != system(c))
1512                             {
1513                                 gmx_fatal(FARGS,"Error executing command: %s",c);
1514                             }
1515 #endif
1516                         }
1517                         outframe++;
1518                     }
1519                 }
1520                 frame++;
1521                 bHaveNextFrame=read_next_frame(oenv,status,&fr);
1522             } while (!(bTDump && bDumpFrame) && bHaveNextFrame);
1523         }
1524
1525         if (!bHaveFirstFrame || (bTDump && !bDumpFrame))
1526             fprintf(stderr,"\nWARNING no output, "
1527                     "last frame read at t=%g\n",fr.time);
1528         fprintf(stderr,"\n");
1529
1530         close_trj(status);
1531         sfree(outf_base);
1532
1533         if (bRmPBC)
1534           gmx_rmpbc_done(gpbc);
1535
1536         if (trxout)
1537             close_trx(trxout);
1538         else if (out != NULL)
1539             ffclose(out);
1540         if (bSubTraj) {
1541             for(i=0; (i<clust->clust->nr); i++)
1542                 if (clust_status[i] )
1543                     close_trx(clust_status[i]);
1544         }
1545     }
1546
1547     do_view(oenv,out_file,NULL);
1548
1549     thanx(stderr);
1550
1551     return 0;
1552 }