src/tools/gmx_trjconv.c

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