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