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