src/tools/gmx_trjconv.c

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