2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
49 #include "gromacs/fileio/gmxfio.h"
50 #include "gromacs/fileio/tpxio.h"
51 #include "gromacs/fileio/trxio.h"
52 #include "gromacs/fileio/trnio.h"
53 #include "gromacs/commandline/pargs.h"
54 #include "gromacs/fileio/futil.h"
55 #include "gromacs/fileio/pdbio.h"
56 #include "gromacs/fileio/confio.h"
60 #include "gromacs/fileio/xtcio.h"
63 #include "gromacs/fileio/g87io.h"
75 euSel, euRect, euTric, euCompact, euNR
79 static void calc_pbc_cluster(int ecenter, int nrefat, t_topology *top, int ePBC,
80 rvec x[], atom_id index[], matrix box)
82 int m, i, j, j0, j1, jj, ai, aj;
85 rvec dx, xtest, box_center;
86 int nmol, imol_center;
88 gmx_bool *bMol, *bTmp;
89 rvec *m_com, *m_shift;
97 calc_box_center(ecenter, box, box_center);
99 /* Initiate the pbc structure */
100 memset(&pbc, 0, sizeof(pbc));
101 set_pbc(&pbc, ePBC, box);
103 /* Convert atom index to molecular */
105 molind = top->mols.index;
111 snew(bTmp, top->atoms.nr);
113 for (i = 0; (i < nrefat); i++)
115 /* Mark all molecules in the index */
118 /* Binary search assuming the molecules are sorted */
123 if (ai < molind[j0+1])
127 else if (ai >= molind[j1])
134 if (ai < molind[jj+1])
146 /* Double check whether all atoms in all molecules that are marked are part
147 * of the cluster. Simultaneously compute the center of geometry.
149 min_dist2 = 10*sqr(trace(box));
152 for (i = 0; i < nmol; i++)
154 for (j = molind[i]; j < molind[i+1]; j++)
156 if (bMol[i] && !bTmp[j])
158 gmx_fatal(FARGS, "Molecule %d marked for clustering but not atom %d in it - check your index!", i+1, j+1);
160 else if (!bMol[i] && bTmp[j])
162 gmx_fatal(FARGS, "Atom %d marked for clustering but not molecule %d - this is an internal error...", j+1, i+1);
166 /* Make molecule whole, move 2nd and higher atom to same periodicity as 1st atom in molecule */
169 pbc_dx(&pbc, x[j], x[j-1], dx);
170 rvec_add(x[j-1], dx, x[j]);
172 /* Compute center of geometry of molecule - m_com[i] was zeroed when we did snew() on it! */
173 rvec_inc(m_com[i], x[j]);
178 /* Normalize center of geometry */
179 fac = 1.0/(molind[i+1]-molind[i]);
180 for (m = 0; (m < DIM); m++)
184 /* Determine which molecule is closest to the center of the box */
185 pbc_dx(&pbc, box_center, m_com[i], dx);
186 tmp_r2 = iprod(dx, dx);
188 if (tmp_r2 < min_dist2)
193 cluster[ncluster++] = i;
200 fprintf(stderr, "No molecules selected in the cluster\n");
203 else if (imol_center == -1)
205 fprintf(stderr, "No central molecules could be found\n");
210 added[nadded++] = imol_center;
211 bMol[imol_center] = FALSE;
213 while (nadded < ncluster)
215 /* Find min distance between cluster molecules and those remaining to be added */
216 min_dist2 = 10*sqr(trace(box));
219 /* Loop over added mols */
220 for (i = 0; i < nadded; i++)
223 /* Loop over all mols */
224 for (j = 0; j < ncluster; j++)
227 /* check those remaining to be added */
230 pbc_dx(&pbc, m_com[aj], m_com[ai], dx);
231 tmp_r2 = iprod(dx, dx);
232 if (tmp_r2 < min_dist2)
242 /* Add the best molecule */
243 added[nadded++] = jmin;
245 /* Calculate the shift from the ai molecule */
246 pbc_dx(&pbc, m_com[jmin], m_com[imin], dx);
247 rvec_add(m_com[imin], dx, xtest);
248 rvec_sub(xtest, m_com[jmin], m_shift[jmin]);
249 rvec_inc(m_com[jmin], m_shift[jmin]);
251 for (j = molind[jmin]; j < molind[jmin+1]; j++)
253 rvec_inc(x[j], m_shift[jmin]);
255 fprintf(stdout, "\rClustering iteration %d of %d...", nadded, ncluster);
265 fprintf(stdout, "\n");
268 static void put_molecule_com_in_box(int unitcell_enum, int ecenter,
270 int natoms, t_atom atom[],
271 int ePBC, matrix box, rvec x[])
275 rvec com, new_com, shift, dx, box_center;
280 calc_box_center(ecenter, box, box_center);
281 set_pbc(&pbc, ePBC, box);
284 gmx_fatal(FARGS, "There are no molecule descriptions. I need a .tpr file for this pbc option.");
286 for (i = 0; (i < mols->nr); i++)
291 for (j = mols->index[i]; (j < mols->index[i+1] && j < natoms); j++)
294 for (d = 0; d < DIM; d++)
300 /* calculate final COM */
301 svmul(1.0/mtot, com, com);
303 /* check if COM is outside box */
304 copy_rvec(com, new_com);
305 switch (unitcell_enum)
308 put_atoms_in_box(ePBC, box, 1, &new_com);
311 put_atoms_in_triclinic_unitcell(ecenter, box, 1, &new_com);
314 put_atoms_in_compact_unitcell(ePBC, ecenter, box, 1, &new_com);
317 rvec_sub(new_com, com, shift);
318 if (norm2(shift) > 0)
322 fprintf(debug, "\nShifting position of molecule %d "
323 "by %8.3f %8.3f %8.3f\n", i+1,
324 shift[XX], shift[YY], shift[ZZ]);
326 for (j = mols->index[i]; (j < mols->index[i+1] && j < natoms); j++)
328 rvec_inc(x[j], shift);
334 static void put_residue_com_in_box(int unitcell_enum, int ecenter,
335 int natoms, t_atom atom[],
336 int ePBC, matrix box, rvec x[])
338 atom_id i, j, res_start, res_end, res_nat;
342 rvec box_center, com, new_com, shift;
344 calc_box_center(ecenter, box, box_center);
350 for (i = 0; i < natoms+1; i++)
352 if (i == natoms || (presnr != atom[i].resind && presnr != NOTSET))
354 /* calculate final COM */
356 res_nat = res_end - res_start;
357 svmul(1.0/mtot, com, com);
359 /* check if COM is outside box */
360 copy_rvec(com, new_com);
361 switch (unitcell_enum)
364 put_atoms_in_box(ePBC, box, 1, &new_com);
367 put_atoms_in_triclinic_unitcell(ecenter, box, 1, &new_com);
370 put_atoms_in_compact_unitcell(ePBC, ecenter, box, 1, &new_com);
373 rvec_sub(new_com, com, shift);
378 fprintf(debug, "\nShifting position of residue %d (atoms %u-%u) "
379 "by %g,%g,%g\n", atom[res_start].resind+1,
380 res_start+1, res_end+1, shift[XX], shift[YY], shift[ZZ]);
382 for (j = res_start; j < res_end; j++)
384 rvec_inc(x[j], shift);
390 /* remember start of new residue */
397 for (d = 0; d < DIM; d++)
403 presnr = atom[i].resind;
408 static void center_x(int ecenter, rvec x[], matrix box, int n, int nc, atom_id ci[])
411 rvec cmin, cmax, box_center, dx;
415 copy_rvec(x[ci[0]], cmin);
416 copy_rvec(x[ci[0]], cmax);
417 for (i = 0; i < nc; i++)
420 for (m = 0; m < DIM; m++)
422 if (x[ai][m] < cmin[m])
426 else if (x[ai][m] > cmax[m])
432 calc_box_center(ecenter, box, box_center);
433 for (m = 0; m < DIM; m++)
435 dx[m] = box_center[m]-(cmin[m]+cmax[m])*0.5;
438 for (i = 0; i < n; i++)
445 static void mk_filenm(char *base, const char *ext, int ndigit, int file_nr,
451 strcpy(out_file, base);
462 strncat(out_file, "00000000000", ndigit-nd);
464 sprintf(nbuf, "%d.", file_nr);
465 strcat(out_file, nbuf);
466 strcat(out_file, ext);
469 void check_trn(const char *fn)
471 if ((fn2ftp(fn) != efTRJ) && (fn2ftp(fn) != efTRR))
473 gmx_fatal(FARGS, "%s is not a trajectory file, exiting\n", fn);
477 #ifndef GMX_NATIVE_WINDOWS
478 void do_trunc(const char *fn, real t0)
491 gmx_fatal(FARGS, "You forgot to set the truncation time");
494 /* Check whether this is a .trj file */
497 in = open_trn(fn, "r");
498 fp = gmx_fio_getfp(in);
501 fprintf(stderr, "Sorry, can not trunc %s, truncation of this filetype is not supported\n", fn);
507 fpos = gmx_fio_ftell(in);
509 while (!bStop && fread_trnheader(in, &sh, &bOK))
511 fread_htrn(in, &sh, NULL, NULL, NULL, NULL);
512 fpos = gmx_ftell(fp);
516 gmx_fseek(fp, fpos, SEEK_SET);
522 fprintf(stderr, "Do you REALLY want to truncate this trajectory (%s) at:\n"
523 "frame %d, time %g, bytes %ld ??? (type YES if so)\n",
524 fn, j, t, (long int)fpos);
525 if (1 != scanf("%s", yesno))
527 gmx_fatal(FARGS, "Error reading user input");
529 if (strcmp(yesno, "YES") == 0)
531 fprintf(stderr, "Once again, I'm gonna DO this...\n");
533 if (0 != truncate(fn, fpos))
535 gmx_fatal(FARGS, "Error truncating file %s", fn);
540 fprintf(stderr, "Ok, I'll forget about it\n");
545 fprintf(stderr, "Already at end of file (t=%g)...\n", t);
552 int gmx_trjconv(int argc, char *argv[])
554 const char *desc[] = {
555 "[THISMODULE] can convert trajectory files in many ways:[BR]",
556 "[BB]1.[bb] from one format to another[BR]",
557 "[BB]2.[bb] select a subset of atoms[BR]",
558 "[BB]3.[bb] change the periodicity representation[BR]",
559 "[BB]4.[bb] keep multimeric molecules together[BR]",
560 "[BB]5.[bb] center atoms in the box[BR]",
561 "[BB]6.[bb] fit atoms to reference structure[BR]",
562 "[BB]7.[bb] reduce the number of frames[BR]",
563 "[BB]8.[bb] change the timestamps of the frames ",
564 "([TT]-t0[tt] and [TT]-timestep[tt])[BR]",
565 "[BB]9.[bb] cut the trajectory in small subtrajectories according",
566 "to information in an index file. This allows subsequent analysis of",
567 "the subtrajectories that could, for example, be the result of a",
568 "cluster analysis. Use option [TT]-sub[tt].",
569 "This assumes that the entries in the index file are frame numbers and",
570 "dumps each group in the index file to a separate trajectory file.[BR]",
571 "[BB]10.[bb] select frames within a certain range of a quantity given",
572 "in an [TT].xvg[tt] file.[PAR]",
574 "[gmx-trjcat] is better suited for concatenating multiple trajectory files.",
577 "Currently seven formats are supported for input and output:",
578 "[TT].xtc[tt], [TT].trr[tt], [TT].trj[tt], [TT].gro[tt], [TT].g96[tt],",
579 "[TT].pdb[tt] and [TT].g87[tt].",
580 "The file formats are detected from the file extension.",
581 "The precision of [TT].xtc[tt] and [TT].gro[tt] output is taken from the",
582 "input file for [TT].xtc[tt], [TT].gro[tt] and [TT].pdb[tt],",
583 "and from the [TT]-ndec[tt] option for other input formats. The precision",
584 "is always taken from [TT]-ndec[tt], when this option is set.",
585 "All other formats have fixed precision. [TT].trr[tt] and [TT].trj[tt]",
586 "output can be single or double precision, depending on the precision",
587 "of the [THISMODULE] binary.",
588 "Note that velocities are only supported in",
589 "[TT].trr[tt], [TT].trj[tt], [TT].gro[tt] and [TT].g96[tt] files.[PAR]",
591 "Option [TT]-app[tt] can be used to",
592 "append output to an existing trajectory file.",
593 "No checks are performed to ensure integrity",
594 "of the resulting combined trajectory file.[PAR]",
596 "Option [TT]-sep[tt] can be used to write every frame to a separate",
597 "[TT].gro, .g96[tt] or [TT].pdb[tt] file. By default, all frames all written to one file.",
598 "[TT].pdb[tt] files with all frames concatenated can be viewed with",
599 "[TT]rasmol -nmrpdb[tt].[PAR]",
601 "It is possible to select part of your trajectory and write it out",
602 "to a new trajectory file in order to save disk space, e.g. for leaving",
603 "out the water from a trajectory of a protein in water.",
604 "[BB]ALWAYS[bb] put the original trajectory on tape!",
605 "We recommend to use the portable [TT].xtc[tt] format for your analysis",
606 "to save disk space and to have portable files.[PAR]",
608 "There are two options for fitting the trajectory to a reference",
609 "either for essential dynamics analysis, etc.",
610 "The first option is just plain fitting to a reference structure",
611 "in the structure file. The second option is a progressive fit",
612 "in which the first timeframe is fitted to the reference structure ",
613 "in the structure file to obtain and each subsequent timeframe is ",
614 "fitted to the previously fitted structure. This way a continuous",
615 "trajectory is generated, which might not be the case when using the",
616 "regular fit method, e.g. when your protein undergoes large",
617 "conformational transitions.[PAR]",
619 "Option [TT]-pbc[tt] sets the type of periodic boundary condition",
621 "[TT]* mol[tt] puts the center of mass of molecules in the box,",
622 "and requires a run input file to be supplied with [TT]-s[tt].[BR]",
623 "[TT]* res[tt] puts the center of mass of residues in the box.[BR]",
624 "[TT]* atom[tt] puts all the atoms in the box.[BR]",
625 "[TT]* nojump[tt] checks if atoms jump across the box and then puts",
626 "them back. This has the effect that all molecules",
627 "will remain whole (provided they were whole in the initial",
628 "conformation). [BB]Note[bb] that this ensures a continuous trajectory but",
629 "molecules may diffuse out of the box. The starting configuration",
630 "for this procedure is taken from the structure file, if one is",
631 "supplied, otherwise it is the first frame.[BR]",
632 "[TT]* cluster[tt] clusters all the atoms in the selected index",
633 "such that they are all closest to the center of mass of the cluster,",
634 "which is iteratively updated. [BB]Note[bb] that this will only give meaningful",
635 "results if you in fact have a cluster. Luckily that can be checked",
636 "afterwards using a trajectory viewer. Note also that if your molecules",
637 "are broken this will not work either.[BR]",
638 "The separate option [TT]-clustercenter[tt] can be used to specify an",
639 "approximate center for the cluster. This is useful e.g. if you have",
640 "two big vesicles, and you want to maintain their relative positions.[BR]",
641 "[TT]* whole[tt] only makes broken molecules whole.[PAR]",
643 "Option [TT]-ur[tt] sets the unit cell representation for options",
644 "[TT]mol[tt], [TT]res[tt] and [TT]atom[tt] of [TT]-pbc[tt].",
645 "All three options give different results for triclinic boxes and",
646 "identical results for rectangular boxes.",
647 "[TT]rect[tt] is the ordinary brick shape.",
648 "[TT]tric[tt] is the triclinic unit cell.",
649 "[TT]compact[tt] puts all atoms at the closest distance from the center",
650 "of the box. This can be useful for visualizing e.g. truncated octahedra",
651 "or rhombic dodecahedra. The center for options [TT]tric[tt] and [TT]compact[tt]",
652 "is [TT]tric[tt] (see below), unless the option [TT]-boxcenter[tt]",
653 "is set differently.[PAR]",
655 "Option [TT]-center[tt] centers the system in the box. The user can",
656 "select the group which is used to determine the geometrical center.",
657 "Option [TT]-boxcenter[tt] sets the location of the center of the box",
658 "for options [TT]-pbc[tt] and [TT]-center[tt]. The center options are:",
659 "[TT]tric[tt]: half of the sum of the box vectors,",
660 "[TT]rect[tt]: half of the box diagonal,",
661 "[TT]zero[tt]: zero.",
662 "Use option [TT]-pbc mol[tt] in addition to [TT]-center[tt] when you",
663 "want all molecules in the box after the centering.[PAR]",
665 "It is not always possible to use combinations of [TT]-pbc[tt],",
666 "[TT]-fit[tt], [TT]-ur[tt] and [TT]-center[tt] to do exactly what",
667 "you want in one call to [THISMODULE]. Consider using multiple",
668 "calls, and check out the GROMACS website for suggestions.[PAR]",
670 "With [TT]-dt[tt], it is possible to reduce the number of ",
671 "frames in the output. This option relies on the accuracy of the times",
672 "in your input trajectory, so if these are inaccurate use the",
673 "[TT]-timestep[tt] option to modify the time (this can be done",
674 "simultaneously). For making smooth movies, the program [gmx-filter]",
675 "can reduce the number of frames while using low-pass frequency",
676 "filtering, this reduces aliasing of high frequency motions.[PAR]",
678 "Using [TT]-trunc[tt] [THISMODULE] can truncate [TT].trj[tt] in place, i.e.",
679 "without copying the file. This is useful when a run has crashed",
680 "during disk I/O (i.e. full disk), or when two contiguous",
681 "trajectories must be concatenated without having double frames.[PAR]",
683 "Option [TT]-dump[tt] can be used to extract a frame at or near",
684 "one specific time from your trajectory.[PAR]",
686 "Option [TT]-drop[tt] reads an [TT].xvg[tt] file with times and values.",
687 "When options [TT]-dropunder[tt] and/or [TT]-dropover[tt] are set,",
688 "frames with a value below and above the value of the respective options",
689 "will not be written."
705 const char *pbc_opt[epNR + 1] =
707 NULL, "none", "mol", "res", "atom", "nojump", "cluster", "whole",
712 const char *unitcell_opt[euNR+1] =
713 { NULL, "rect", "tric", "compact", NULL };
717 ecSel, ecTric, ecRect, ecZero, ecNR
719 const char *center_opt[ecNR+1] =
720 { NULL, "tric", "rect", "zero", NULL };
726 efSel, efNone, efFit, efFitXY, efReset, efResetXY, efPFit, efNR
728 const char *fit[efNR + 1] =
730 NULL, "none", "rot+trans", "rotxy+transxy", "translation", "transxy",
734 static gmx_bool bAppend = FALSE, bSeparate = FALSE, bVels = TRUE, bForce = FALSE, bCONECT = FALSE;
735 static gmx_bool bCenter = FALSE;
736 static int skip_nr = 1, ndec = 3, nzero = 0;
737 static real tzero = 0, delta_t = 0, timestep = 0, ttrunc = -1, tdump = -1, split_t = 0;
738 static rvec newbox = {0, 0, 0}, shift = {0, 0, 0}, trans = {0, 0, 0};
739 static char *exec_command = NULL;
740 static real dropunder = 0, dropover = 0;
741 static gmx_bool bRound = FALSE;
746 { "-skip", FALSE, etINT,
747 { &skip_nr }, "Only write every nr-th frame" },
748 { "-dt", FALSE, etTIME,
750 "Only write frame when t MOD dt = first time (%t)" },
751 { "-round", FALSE, etBOOL,
752 { &bRound }, "Round measurements to nearest picosecond"},
753 { "-dump", FALSE, etTIME,
754 { &tdump }, "Dump frame nearest specified time (%t)" },
755 { "-t0", FALSE, etTIME,
757 "Starting time (%t) (default: don't change)" },
758 { "-timestep", FALSE, etTIME,
760 "Change time step between input frames (%t)" },
761 { "-pbc", FALSE, etENUM,
763 "PBC treatment (see help text for full description)" },
764 { "-ur", FALSE, etENUM,
765 { unitcell_opt }, "Unit-cell representation" },
766 { "-center", FALSE, etBOOL,
767 { &bCenter }, "Center atoms in box" },
768 { "-boxcenter", FALSE, etENUM,
769 { center_opt }, "Center for -pbc and -center" },
770 { "-box", FALSE, etRVEC,
772 "Size for new cubic box (default: read from input)" },
773 { "-trans", FALSE, etRVEC,
775 "All coordinates will be translated by trans. This "
776 "can advantageously be combined with -pbc mol -ur "
778 { "-shift", FALSE, etRVEC,
780 "All coordinates will be shifted by framenr*shift" },
781 { "-fit", FALSE, etENUM,
783 "Fit molecule to ref structure in the structure file" },
784 { "-ndec", FALSE, etINT,
786 "Precision for .xtc and .gro writing in number of "
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,
795 "Truncate input trajectory file after this time (%t)" },
797 { "-exec", FALSE, etSTR,
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,
805 "Start writing new file when t MOD split = first "
807 { "-sep", FALSE, etBOOL,
809 "Write each frame to a separate .gro, .g96 or .pdb "
811 { "-nzero", FALSE, etINT,
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,
821 "Add conect records when writing [TT].pdb[tt] files. Useful "
822 "for visualization of non-standard molecules, e.g. "
823 "coarse grained ones" }
825 #define NPA asize(pa)
828 t_trxstatus *trxout = NULL;
830 int ftp, ftpin = 0, file_nr;
831 t_trxframe fr, frout;
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;
839 gmx_conect gc = NULL;
841 t_atoms *atoms = NULL, useatoms;
843 atom_id *index, *cindex;
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;
854 int *nfwritten = NULL;
855 int ndrop = 0, ncol, drop0 = 0, drop1 = 0, dropuse = 0;
857 real tshift = 0, t0 = -1, dt = 0.001, prec;
858 gmx_bool bFit, bFitXY, bPFit, bReset;
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];
873 gmx_bool bWarnCompact = FALSE;
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 }
886 #define NFILE asize(fnm)
888 if (!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,
897 top_file = ftp2fn(efTPS, NFILE, fnm);
900 /* Check command line */
901 in_file = opt2fn("-f", NFILE, fnm);
905 #ifndef GMX_NATIVE_WINDOWS
906 do_trunc(in_file, ttrunc);
911 /* mark active cmdline options */
912 bSetBox = opt2parg_bSet("-box", NPA, pa);
913 bSetTime = opt2parg_bSet("-t0", NPA, pa);
914 bSetPrec = opt2parg_bSet("-ndec", NPA, pa);
915 bSetUR = opt2parg_bSet("-ur", NPA, pa);
916 bExec = opt2parg_bSet("-exec", NPA, pa);
917 bTimeStep = opt2parg_bSet("-timestep", NPA, pa);
918 bTDump = opt2parg_bSet("-dump", NPA, pa);
919 bDropUnder = opt2parg_bSet("-dropunder", NPA, pa);
920 bDropOver = opt2parg_bSet("-dropover", NPA, pa);
921 bTrans = opt2parg_bSet("-trans", NPA, pa);
922 bSplit = (split_t != 0);
924 /* parse enum options */
925 fit_enum = nenum(fit);
926 bFit = (fit_enum == efFit || fit_enum == efFitXY);
927 bFitXY = fit_enum == efFitXY;
928 bReset = (fit_enum == efReset || fit_enum == efResetXY);
929 bPFit = fit_enum == efPFit;
930 pbc_enum = nenum(pbc_opt);
931 bPBCWhole = pbc_enum == epWhole;
932 bPBCcomRes = pbc_enum == epComRes;
933 bPBCcomMol = pbc_enum == epComMol;
934 bPBCcomAtom = pbc_enum == epComAtom;
935 bNoJump = pbc_enum == epNojump;
936 bCluster = pbc_enum == epCluster;
937 bPBC = pbc_enum != epNone;
938 unitcell_enum = nenum(unitcell_opt);
939 ecenter = nenum(center_opt) - ecTric;
941 /* set and check option dependencies */
944 bFit = TRUE; /* for pfit, fit *must* be set */
948 bReset = TRUE; /* for fit, reset *must* be set */
953 nfitdim = (fit_enum == efFitXY || fit_enum == efResetXY) ? 2 : 3;
955 bRmPBC = bFit || bPBCWhole || bPBCcomRes || bPBCcomMol;
959 if (!(bPBCcomRes || bPBCcomMol || bPBCcomAtom))
962 "WARNING: Option for unitcell representation (-ur %s)\n"
963 " only has effect in combination with -pbc %s, %s or %s.\n"
964 " Ingoring unitcell representation.\n\n",
965 unitcell_opt[0], pbc_opt[2], pbc_opt[3], pbc_opt[4]);
971 gmx_fatal(FARGS, "PBC condition treatment does not work together with rotational fit.\n"
972 "Please do the PBC condition treatment first and then run trjconv in a second step\n"
973 "for the rotational fit.\n"
974 "First doing the rotational fit and then doing the PBC treatment gives incorrect\n"
978 /* ndec is in nr of decimal places, prec is a multiplication factor: */
980 for (i = 0; i < ndec; i++)
985 bIndex = ftp2bSet(efNDX, NFILE, fnm);
988 /* Determine output type */
989 out_file = opt2fn("-o", NFILE, fnm);
990 ftp = fn2ftp(out_file);
991 fprintf(stderr, "Will write %s: %s\n", ftp2ext(ftp), ftp2desc(ftp));
992 bNeedPrec = (ftp == efXTC || ftp == efGRO);
995 /* check if velocities are possible in input and output files */
996 ftpin = fn2ftp(in_file);
997 bVels = (ftp == efTRR || ftp == efTRJ || ftp == efGRO || ftp == efG96)
998 && (ftpin == efTRR || ftpin == efTRJ || ftpin == efGRO || ftpin == efG96 ||
1001 if (bSeparate || bSplit)
1003 outf_ext = strrchr(out_file, '.');
1004 if (outf_ext == NULL)
1006 gmx_fatal(FARGS, "Output file name '%s' does not contain a '.'", out_file);
1008 outf_base = strdup(out_file);
1009 outf_base[outf_ext - out_file] = '\0';
1012 bSubTraj = opt2bSet("-sub", NFILE, fnm);
1015 if ((ftp != efXTC) && (ftp != efTRR))
1017 /* It seems likely that other trajectory file types
1018 * could work here. */
1019 gmx_fatal(FARGS, "Can only use the sub option with output file types "
1022 clust = cluster_index(NULL, opt2fn("-sub", NFILE, fnm));
1024 /* Check for number of files disabled, as FOPEN_MAX is not the correct
1025 * number to check for. In my linux box it is only 16.
1027 if (0 && (clust->clust->nr > FOPEN_MAX-4))
1029 gmx_fatal(FARGS, "Can not open enough (%d) files to write all the"
1030 " trajectories.\ntry splitting the index file in %d parts.\n"
1032 clust->clust->nr, 1+clust->clust->nr/FOPEN_MAX, FOPEN_MAX);
1034 gmx_warning("The -sub option could require as many open output files as there are\n"
1035 "index groups in the file (%d). If you get I/O errors opening new files,\n"
1036 "try reducing the number of index groups in the file, and perhaps\n"
1037 "using trjconv -sub several times on different chunks of your index file.\n",
1040 snew(clust_status, clust->clust->nr);
1041 snew(clust_status_id, clust->clust->nr);
1042 snew(nfwritten, clust->clust->nr);
1043 for (i = 0; (i < clust->clust->nr); i++)
1045 clust_status[i] = NULL;
1046 clust_status_id[i] = -1;
1048 bSeparate = bSplit = FALSE;
1055 /* Determine whether to read a topology */
1056 bTPS = (ftp2bSet(efTPS, NFILE, fnm) ||
1057 bRmPBC || bReset || bPBCcomMol || bCluster ||
1058 (ftp == efGRO) || (ftp == efPDB) || bCONECT);
1060 /* Determine if when can read index groups */
1061 bIndex = (bIndex || bTPS);
1065 read_tps_conf(top_file, top_title, &top, &ePBC, &xp, NULL, top_box,
1066 bReset || bPBCcomRes);
1069 if (0 == top.mols.nr && (bCluster || bPBCcomMol))
1071 gmx_fatal(FARGS, "Option -pbc %s requires a .tpr file for the -s option", pbc_opt[pbc_enum]);
1074 /* top_title is only used for gro and pdb,
1075 * the header in such a file is top_title t= ...
1076 * to prevent a double t=, remove it from top_title
1078 if ((charpt = strstr(top_title, " t= ")))
1085 gc = gmx_conect_generate(&top);
1089 gpbc = gmx_rmpbc_init(&top.idef, ePBC, top.atoms.nr);
1093 /* get frame number index */
1095 if (opt2bSet("-fr", NFILE, fnm))
1097 printf("Select groups of frame number indices:\n");
1098 rd_index(opt2fn("-fr", NFILE, fnm), 1, &nrfri, (atom_id **)&frindex, &frname);
1101 for (i = 0; i < nrfri; i++)
1103 fprintf(debug, "frindex[%4d]=%4d\n", i, frindex[i]);
1108 /* get index groups etc. */
1111 printf("Select group for %s fit\n",
1112 bFit ? "least squares" : "translational");
1113 get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm),
1114 1, &ifit, &ind_fit, &gn_fit);
1120 gmx_fatal(FARGS, "Need at least 2 atoms to fit!\n");
1124 fprintf(stderr, "WARNING: fitting with only 2 atoms is not unique\n");
1130 printf("Select group for clustering\n");
1131 get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm),
1132 1, &ifit, &ind_fit, &gn_fit);
1139 printf("Select group for centering\n");
1140 get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm),
1141 1, &ncent, &cindex, &grpnm);
1143 printf("Select group for output\n");
1144 get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm),
1145 1, &nout, &index, &grpnm);
1149 /* no index file, so read natoms from TRX */
1150 if (!read_first_frame(oenv, &status, in_file, &fr, TRX_DONT_SKIP))
1152 gmx_fatal(FARGS, "Could not read a frame from %s", in_file);
1157 snew(index, natoms);
1158 for (i = 0; i < natoms; i++)
1172 snew(w_rls, atoms->nr);
1173 for (i = 0; (i < ifit); i++)
1175 w_rls[ind_fit[i]] = atoms->atom[ind_fit[i]].m;
1178 /* Restore reference structure and set to origin,
1179 store original location (to put structure back) */
1182 gmx_rmpbc(gpbc, top.atoms.nr, top_box, xp);
1184 copy_rvec(xp[index[0]], x_shift);
1185 reset_x_ndim(nfitdim, ifit, ind_fit, atoms->nr, NULL, xp, w_rls);
1186 rvec_dec(x_shift, xp[index[0]]);
1190 clear_rvec(x_shift);
1193 if (bDropUnder || bDropOver)
1195 /* Read the .xvg file with the drop values */
1196 fprintf(stderr, "\nReading drop file ...");
1197 ndrop = read_xvg(opt2fn("-drop", NFILE, fnm), &dropval, &ncol);
1198 fprintf(stderr, " %d time points\n", ndrop);
1199 if (ndrop == 0 || ncol < 2)
1201 gmx_fatal(FARGS, "Found no data points in %s",
1202 opt2fn("-drop", NFILE, fnm));
1208 /* Make atoms struct for output in GRO or PDB files */
1209 if ((ftp == efGRO) || ((ftp == efG96) && bTPS) || (ftp == efPDB))
1211 /* get memory for stuff to go in .pdb file */
1212 init_t_atoms(&useatoms, atoms->nr, FALSE);
1213 sfree(useatoms.resinfo);
1214 useatoms.resinfo = atoms->resinfo;
1215 for (i = 0; (i < nout); i++)
1217 useatoms.atomname[i] = atoms->atomname[index[i]];
1218 useatoms.atom[i] = atoms->atom[index[i]];
1219 useatoms.nres = max(useatoms.nres, useatoms.atom[i].resind+1);
1223 /* select what to read */
1224 if (ftp == efTRR || ftp == efTRJ)
1234 flags = flags | TRX_READ_V;
1238 flags = flags | TRX_READ_F;
1241 /* open trx file for reading */
1242 bHaveFirstFrame = read_first_frame(oenv, &status, in_file, &fr, flags);
1245 fprintf(stderr, "\nPrecision of %s is %g (nm)\n", in_file, 1/fr.prec);
1249 if (bSetPrec || !fr.bPrec)
1251 fprintf(stderr, "\nSetting output precision to %g (nm)\n", 1/prec);
1255 fprintf(stderr, "Using output precision of %g (nm)\n", 1/prec);
1259 if (bHaveFirstFrame)
1261 set_trxframe_ePBC(&fr, ePBC);
1267 tshift = tzero-fr.time;
1274 /* open output for writing */
1275 if ((bAppend) && (gmx_fexist(out_file)))
1277 strcpy(filemode, "a");
1278 fprintf(stderr, "APPENDING to existing file %s\n", out_file);
1282 strcpy(filemode, "w");
1291 if (!bSplit && !bSubTraj)
1293 trxout = open_trx(out_file, filemode);
1299 if (( !bSeparate && !bSplit ) && !bSubTraj)
1301 out = ffopen(out_file, filemode);
1309 /* check if index is meaningful */
1310 for (i = 0; i < nout; i++)
1312 if (index[i] >= natoms)
1315 "Index[%d] %d is larger than the number of atoms in the\n"
1316 "trajectory file (%d). There is a mismatch in the contents\n"
1317 "of your -f, -s and/or -n files.", i, index[i]+1, natoms);
1319 bCopy = bCopy || (i != index[i]);
1337 fprintf(gmx_fio_getfp(trx_get_fileio(trxout)),
1338 "Generated by %s. #atoms=%d, a BOX is"
1339 " stored in this file.\n", ShortProgram(), nout);
1342 /* Start the big loop over frames */
1349 /* Main loop over frames */
1360 /*if (frame >= clust->clust->nra)
1361 gmx_fatal(FARGS,"There are more frames in the trajectory than in the cluster index file\n");*/
1362 if (frame > clust->maxframe)
1368 my_clust = clust->inv_clust[frame];
1370 if ((my_clust < 0) || (my_clust >= clust->clust->nr) ||
1371 (my_clust == NO_ATID))
1379 /* generate new box */
1381 for (m = 0; m < DIM; m++)
1383 fr.box[m][m] = newbox[m];
1389 for (i = 0; i < natoms; i++)
1391 rvec_inc(fr.x[i], trans);
1397 /* determine timestep */
1410 /* This is not very elegant, as one can not dump a frame after
1411 * a timestep with is more than twice as small as the first one. */
1412 bDumpFrame = (fr.time > tdump-0.5*dt) && (fr.time <= tdump+0.5*dt);
1419 /* determine if an atom jumped across the box and reset it if so */
1420 if (bNoJump && (bTPS || frame != 0))
1422 for (d = 0; d < DIM; d++)
1424 hbox[d] = 0.5*fr.box[d][d];
1426 for (i = 0; i < natoms; i++)
1430 rvec_dec(fr.x[i], x_shift);
1432 for (m = DIM-1; m >= 0; m--)
1436 while (fr.x[i][m]-xp[i][m] <= -hbox[m])
1438 for (d = 0; d <= m; d++)
1440 fr.x[i][d] += fr.box[m][d];
1443 while (fr.x[i][m]-xp[i][m] > hbox[m])
1445 for (d = 0; d <= m; d++)
1447 fr.x[i][d] -= fr.box[m][d];
1456 calc_pbc_cluster(ecenter, ifit, &top, ePBC, fr.x, ind_fit, fr.box);
1461 /* Now modify the coords according to the flags,
1462 for normal fit, this is only done for output frames */
1465 gmx_rmpbc_trxfr(gpbc, &fr);
1468 reset_x_ndim(nfitdim, ifit, ind_fit, natoms, NULL, fr.x, w_rls);
1469 do_fit(natoms, w_rls, xp, fr.x);
1472 /* store this set of coordinates for future use */
1473 if (bPFit || bNoJump)
1479 for (i = 0; (i < natoms); i++)
1481 copy_rvec(fr.x[i], xp[i]);
1482 rvec_inc(fr.x[i], x_shift);
1488 /* see if we have a frame from the frame index group */
1489 for (i = 0; i < nrfri && !bDumpFrame; i++)
1491 bDumpFrame = frame == frindex[i];
1494 if (debug && bDumpFrame)
1496 fprintf(debug, "dumping %d\n", frame);
1500 ( ( !bTDump && !frindex && frame % skip_nr == 0 ) || bDumpFrame );
1502 if (bWriteFrame && (bDropUnder || bDropOver))
1504 while (dropval[0][drop1] < fr.time && drop1+1 < ndrop)
1509 if (fabs(dropval[0][drop0] - fr.time)
1510 < fabs(dropval[0][drop1] - fr.time))
1518 if ((bDropUnder && dropval[1][dropuse] < dropunder) ||
1519 (bDropOver && dropval[1][dropuse] > dropover))
1521 bWriteFrame = FALSE;
1531 fr.time = tzero+frame*timestep;
1541 fprintf(stderr, "\nDumping frame at t= %g %s\n",
1542 output_env_conv_time(oenv, fr.time), output_env_get_time_unit(oenv));
1545 /* check for writing at each delta_t */
1546 bDoIt = (delta_t == 0);
1551 bDoIt = bRmod(fr.time, tzero, delta_t);
1555 /* round() is not C89 compatible, so we do this: */
1556 bDoIt = bRmod(floor(fr.time+0.5), floor(tzero+0.5),
1557 floor(delta_t+0.5));
1561 if (bDoIt || bTDump)
1563 /* print sometimes */
1564 if ( ((outframe % SKIP) == 0) || (outframe < SKIP) )
1566 fprintf(stderr, " -> frame %6d time %8.3f \r",
1567 outframe, output_env_conv_time(oenv, fr.time));
1572 /* Now modify the coords according to the flags,
1573 for PFit we did this already! */
1577 gmx_rmpbc_trxfr(gpbc, &fr);
1582 reset_x_ndim(nfitdim, ifit, ind_fit, natoms, NULL, fr.x, w_rls);
1585 do_fit_ndim(nfitdim, natoms, w_rls, xp, fr.x);
1589 for (i = 0; i < natoms; i++)
1591 rvec_inc(fr.x[i], x_shift);
1598 center_x(ecenter, fr.x, fr.box, natoms, ncent, cindex);
1604 switch (unitcell_enum)
1607 put_atoms_in_box(ePBC, fr.box, natoms, fr.x);
1610 put_atoms_in_triclinic_unitcell(ecenter, fr.box, natoms, fr.x);
1613 warn = put_atoms_in_compact_unitcell(ePBC, ecenter, fr.box,
1615 if (warn && !bWarnCompact)
1617 fprintf(stderr, "\n%s\n", warn);
1618 bWarnCompact = TRUE;
1625 put_residue_com_in_box(unitcell_enum, ecenter,
1626 natoms, atoms->atom, ePBC, fr.box, fr.x);
1630 put_molecule_com_in_box(unitcell_enum, ecenter,
1632 natoms, atoms->atom, ePBC, fr.box, fr.x);
1634 /* Copy the input trxframe struct to the output trxframe struct */
1636 frout.bV = (frout.bV && bVels);
1637 frout.bF = (frout.bF && bForce);
1638 frout.natoms = nout;
1639 if (bNeedPrec && (bSetPrec || !fr.bPrec))
1655 for (i = 0; i < nout; i++)
1657 copy_rvec(fr.x[index[i]], frout.x[i]);
1660 copy_rvec(fr.v[index[i]], frout.v[i]);
1664 copy_rvec(fr.f[index[i]], frout.f[i]);
1669 if (opt2parg_bSet("-shift", NPA, pa))
1671 for (i = 0; i < nout; i++)
1673 for (d = 0; d < DIM; d++)
1675 frout.x[i][d] += outframe*shift[d];
1682 bSplitHere = bSplit && bRmod(fr.time, tzero, split_t);
1686 /* round() is not C89 compatible, so we do this: */
1687 bSplitHere = bSplit && bRmod(floor(fr.time+0.5),
1689 floor(split_t+0.5));
1691 if (bSeparate || bSplitHere)
1693 mk_filenm(outf_base, ftp2ext(ftp), nzero, file_nr, out_file2);
1708 trxout = open_trx(out_file2, filemode);
1715 if (clust_status_id[my_clust] == -1)
1717 sprintf(buf, "%s.%s", clust->grpname[my_clust], ftp2ext(ftp));
1718 clust_status[my_clust] = open_trx(buf, "w");
1719 clust_status_id[my_clust] = 1;
1722 else if (clust_status_id[my_clust] == -2)
1724 gmx_fatal(FARGS, "File %s.xtc should still be open (%d open .xtc files)\n" "in order to write frame %d. my_clust = %d",
1725 clust->grpname[my_clust], ntrxopen, frame,
1728 write_trxframe(clust_status[my_clust], &frout, gc);
1729 nfwritten[my_clust]++;
1730 if (nfwritten[my_clust] ==
1731 (clust->clust->index[my_clust+1]-
1732 clust->clust->index[my_clust]))
1734 close_trx(clust_status[my_clust]);
1735 clust_status[my_clust] = NULL;
1736 clust_status_id[my_clust] = -2;
1740 gmx_fatal(FARGS, "Less than zero open .xtc files!");
1747 write_trxframe(trxout, &frout, gc);
1753 sprintf(title, "Generated by trjconv : %s t= %9.5f",
1754 top_title, fr.time);
1755 if (bSeparate || bSplitHere)
1757 out = ffopen(out_file2, "w");
1762 write_hconf_p(out, title, &useatoms, prec2ndec(frout.prec),
1763 frout.x, frout.bV ? frout.v : NULL, frout.box);
1766 fprintf(out, "REMARK GENERATED BY TRJCONV\n");
1767 sprintf(title, "%s t= %9.5f", top_title, frout.time);
1768 /* if reading from pdb, we want to keep the original
1769 model numbering else we write the output frame
1770 number plus one, because model 0 is not allowed in pdb */
1771 if (ftpin == efPDB && fr.bStep && fr.step > model_nr)
1779 write_pdbfile(out, title, &useatoms, frout.x,
1780 frout.ePBC, frout.box, ' ', model_nr, gc, TRUE);
1783 frout.title = title;
1784 if (bSeparate || bTDump)
1786 frout.bTitle = TRUE;
1789 frout.bAtoms = TRUE;
1791 frout.atoms = &useatoms;
1792 frout.bStep = FALSE;
1793 frout.bTime = FALSE;
1797 frout.bTitle = (outframe == 0);
1798 frout.bAtoms = FALSE;
1802 write_g96_conf(out, &frout, -1, NULL);
1811 gmx_fatal(FARGS, "DHE, ftp=%d\n", ftp);
1813 if (bSeparate || bSplitHere)
1818 /* execute command */
1822 sprintf(c, "%s %d", exec_command, file_nr-1);
1823 /*fprintf(stderr,"Executing '%s'\n",c);*/
1824 #ifdef GMX_NO_SYSTEM
1825 printf("Warning-- No calls to system(3) supported on this platform.");
1826 printf("Warning-- Skipping execution of 'system(\"%s\")'.", c);
1830 gmx_fatal(FARGS, "Error executing command: %s", c);
1838 bHaveNextFrame = read_next_frame(oenv, status, &fr);
1840 while (!(bTDump && bDumpFrame) && bHaveNextFrame);
1843 if (!bHaveFirstFrame || (bTDump && !bDumpFrame))
1845 fprintf(stderr, "\nWARNING no output, "
1846 "last frame read at t=%g\n", fr.time);
1848 fprintf(stderr, "\n");
1855 gmx_rmpbc_done(gpbc);
1862 else if (out != NULL)
1868 for (i = 0; (i < clust->clust->nr); i++)
1870 if (clust_status_id[i] >= 0)
1872 close_trx(clust_status[i]);
1878 do_view(oenv, out_file, NULL);