src/gromacs/gmxana/gmx_trjconv.c

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