3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
70 euSel, euRect, euTric, euCompact, euNR
74 static void calc_pbc_cluster(int ecenter, int nrefat, t_topology *top, int ePBC,
75 rvec x[], atom_id index[], matrix box)
77 int m, i, j, j0, j1, jj, ai, aj;
80 rvec dx, xtest, box_center;
81 int nmol, imol_center;
83 gmx_bool *bMol, *bTmp;
84 rvec *m_com, *m_shift;
92 calc_box_center(ecenter, box, box_center);
94 /* Initiate the pbc structure */
95 memset(&pbc, 0, sizeof(pbc));
96 set_pbc(&pbc, ePBC, box);
98 /* Convert atom index to molecular */
100 molind = top->mols.index;
106 snew(bTmp, top->atoms.nr);
108 for (i = 0; (i < nrefat); i++)
110 /* Mark all molecules in the index */
113 /* Binary search assuming the molecules are sorted */
118 if (ai < molind[j0+1])
122 else if (ai >= molind[j1])
129 if (ai < molind[jj+1])
141 /* Double check whether all atoms in all molecules that are marked are part
142 * of the cluster. Simultaneously compute the center of geometry.
144 min_dist2 = 10*sqr(trace(box));
147 for (i = 0; i < nmol; i++)
149 for (j = molind[i]; j < molind[i+1]; j++)
151 if (bMol[i] && !bTmp[j])
153 gmx_fatal(FARGS, "Molecule %d marked for clustering but not atom %d in it - check your index!", i+1, j+1);
155 else if (!bMol[i] && bTmp[j])
157 gmx_fatal(FARGS, "Atom %d marked for clustering but not molecule %d - this is an internal error...", j+1, i+1);
161 /* Make molecule whole, move 2nd and higher atom to same periodicity as 1st atom in molecule */
164 pbc_dx(&pbc, x[j], x[j-1], dx);
165 rvec_add(x[j-1], dx, x[j]);
167 /* Compute center of geometry of molecule - m_com[i] was zeroed when we did snew() on it! */
168 rvec_inc(m_com[i], x[j]);
173 /* Normalize center of geometry */
174 fac = 1.0/(molind[i+1]-molind[i]);
175 for (m = 0; (m < DIM); m++)
179 /* Determine which molecule is closest to the center of the box */
180 pbc_dx(&pbc, box_center, m_com[i], dx);
181 tmp_r2 = iprod(dx, dx);
183 if (tmp_r2 < min_dist2)
188 cluster[ncluster++] = i;
195 fprintf(stderr, "No molecules selected in the cluster\n");
198 else if (imol_center == -1)
200 fprintf(stderr, "No central molecules could be found\n");
205 added[nadded++] = imol_center;
206 bMol[imol_center] = FALSE;
208 while (nadded < ncluster)
210 /* Find min distance between cluster molecules and those remaining to be added */
211 min_dist2 = 10*sqr(trace(box));
214 /* Loop over added mols */
215 for (i = 0; i < nadded; i++)
218 /* Loop over all mols */
219 for (j = 0; j < ncluster; j++)
222 /* check those remaining to be added */
225 pbc_dx(&pbc, m_com[aj], m_com[ai], dx);
226 tmp_r2 = iprod(dx, dx);
227 if (tmp_r2 < min_dist2)
237 /* Add the best molecule */
238 added[nadded++] = jmin;
240 /* Calculate the shift from the ai molecule */
241 pbc_dx(&pbc, m_com[jmin], m_com[imin], dx);
242 rvec_add(m_com[imin], dx, xtest);
243 rvec_sub(xtest, m_com[jmin], m_shift[jmin]);
244 rvec_inc(m_com[jmin], m_shift[jmin]);
246 for (j = molind[jmin]; j < molind[jmin+1]; j++)
248 rvec_inc(x[j], m_shift[jmin]);
250 fprintf(stdout, "\rClustering iteration %d of %d...", nadded, ncluster);
260 fprintf(stdout, "\n");
263 static void put_molecule_com_in_box(int unitcell_enum, int ecenter,
265 int natoms, t_atom atom[],
266 int ePBC, matrix box, rvec x[])
270 rvec com, new_com, shift, dx, box_center;
275 calc_box_center(ecenter, box, box_center);
276 set_pbc(&pbc, ePBC, box);
279 gmx_fatal(FARGS, "There are no molecule descriptions. I need a .tpr file for this pbc option.");
281 for (i = 0; (i < mols->nr); i++)
286 for (j = mols->index[i]; (j < mols->index[i+1] && j < natoms); j++)
289 for (d = 0; d < DIM; d++)
295 /* calculate final COM */
296 svmul(1.0/mtot, com, com);
298 /* check if COM is outside box */
299 copy_rvec(com, new_com);
300 switch (unitcell_enum)
303 put_atoms_in_box(ePBC, box, 1, &new_com);
306 put_atoms_in_triclinic_unitcell(ecenter, box, 1, &new_com);
309 put_atoms_in_compact_unitcell(ePBC, ecenter, box, 1, &new_com);
312 rvec_sub(new_com, com, shift);
313 if (norm2(shift) > 0)
317 fprintf (debug, "\nShifting position of molecule %d "
318 "by %8.3f %8.3f %8.3f\n", i+1, PR_VEC(shift));
320 for (j = mols->index[i]; (j < mols->index[i+1] && j < natoms); j++)
322 rvec_inc(x[j], shift);
328 static void put_residue_com_in_box(int unitcell_enum, int ecenter,
329 int natoms, t_atom atom[],
330 int ePBC, matrix box, rvec x[])
332 atom_id i, j, res_start, res_end, res_nat;
336 rvec box_center, com, new_com, shift;
338 calc_box_center(ecenter, box, box_center);
344 for (i = 0; i < natoms+1; i++)
346 if (i == natoms || (presnr != atom[i].resind && presnr != NOTSET))
348 /* calculate final COM */
350 res_nat = res_end - res_start;
351 svmul(1.0/mtot, com, com);
353 /* check if COM is outside box */
354 copy_rvec(com, new_com);
355 switch (unitcell_enum)
358 put_atoms_in_box(ePBC, box, 1, &new_com);
361 put_atoms_in_triclinic_unitcell(ecenter, box, 1, &new_com);
364 put_atoms_in_compact_unitcell(ePBC, ecenter, box, 1, &new_com);
367 rvec_sub(new_com, com, shift);
372 fprintf (debug, "\nShifting position of residue %d (atoms %u-%u) "
373 "by %g,%g,%g\n", atom[res_start].resind+1,
374 res_start+1, res_end+1, PR_VEC(shift));
376 for (j = res_start; j < res_end; j++)
378 rvec_inc(x[j], shift);
384 /* remember start of new residue */
391 for (d = 0; d < DIM; d++)
397 presnr = atom[i].resind;
402 static void center_x(int ecenter, rvec x[], matrix box, int n, int nc, atom_id ci[])
405 rvec cmin, cmax, box_center, dx;
409 copy_rvec(x[ci[0]], cmin);
410 copy_rvec(x[ci[0]], cmax);
411 for (i = 0; i < nc; i++)
414 for (m = 0; m < DIM; m++)
416 if (x[ai][m] < cmin[m])
420 else if (x[ai][m] > cmax[m])
426 calc_box_center(ecenter, box, box_center);
427 for (m = 0; m < DIM; m++)
429 dx[m] = box_center[m]-(cmin[m]+cmax[m])*0.5;
432 for (i = 0; i < n; i++)
439 static void mk_filenm(char *base, const char *ext, int ndigit, int file_nr,
445 strcpy(out_file, base);
456 strncat(out_file, "00000000000", ndigit-nd);
458 sprintf(nbuf, "%d.", file_nr);
459 strcat(out_file, nbuf);
460 strcat(out_file, ext);
463 void check_trn(const char *fn)
465 if ((fn2ftp(fn) != efTRJ) && (fn2ftp(fn) != efTRR))
467 gmx_fatal(FARGS, "%s is not a trajectory file, exiting\n", fn);
471 #ifndef GMX_NATIVE_WINDOWS
472 void do_trunc(const char *fn, real t0)
485 gmx_fatal(FARGS, "You forgot to set the truncation time");
488 /* Check whether this is a .trj file */
491 in = open_trn(fn, "r");
492 fp = gmx_fio_getfp(in);
495 fprintf(stderr, "Sorry, can not trunc %s, truncation of this filetype is not supported\n", fn);
501 fpos = gmx_fio_ftell(in);
503 while (!bStop && fread_trnheader(in, &sh, &bOK))
505 fread_htrn(in, &sh, NULL, NULL, NULL, NULL);
506 fpos = gmx_ftell(fp);
510 gmx_fseek(fp, fpos, SEEK_SET);
516 fprintf(stderr, "Do you REALLY want to truncate this trajectory (%s) at:\n"
517 "frame %d, time %g, bytes %ld ??? (type YES if so)\n",
518 fn, j, t, (long int)fpos);
519 if (1 != scanf("%s", yesno))
521 gmx_fatal(FARGS, "Error reading user input");
523 if (strcmp(yesno, "YES") == 0)
525 fprintf(stderr, "Once again, I'm gonna DO this...\n");
527 if (0 != truncate(fn, fpos))
529 gmx_fatal(FARGS, "Error truncating file %s", fn);
534 fprintf(stderr, "Ok, I'll forget about it\n");
539 fprintf(stderr, "Already at end of file (t=%g)...\n", t);
546 int gmx_trjconv(int argc, char *argv[])
548 const char *desc[] = {
549 "[TT]trjconv[tt] can convert trajectory files in many ways:[BR]",
550 "[BB]1.[bb] from one format to another[BR]",
551 "[BB]2.[bb] select a subset of atoms[BR]",
552 "[BB]3.[bb] change the periodicity representation[BR]",
553 "[BB]4.[bb] keep multimeric molecules together[BR]",
554 "[BB]5.[bb] center atoms in the box[BR]",
555 "[BB]6.[bb] fit atoms to reference structure[BR]",
556 "[BB]7.[bb] reduce the number of frames[BR]",
557 "[BB]8.[bb] change the timestamps of the frames ",
558 "([TT]-t0[tt] and [TT]-timestep[tt])[BR]",
559 "[BB]9.[bb] cut the trajectory in small subtrajectories according",
560 "to information in an index file. This allows subsequent analysis of",
561 "the subtrajectories that could, for example, be the result of a",
562 "cluster analysis. Use option [TT]-sub[tt].",
563 "This assumes that the entries in the index file are frame numbers and",
564 "dumps each group in the index file to a separate trajectory file.[BR]",
565 "[BB]10.[bb] select frames within a certain range of a quantity given",
566 "in an [TT].xvg[tt] file.[PAR]",
568 "The program [TT]trjcat[tt] is better suited for concatenating multiple trajectory files.",
571 "Currently seven formats are supported for input and output:",
572 "[TT].xtc[tt], [TT].trr[tt], [TT].trj[tt], [TT].gro[tt], [TT].g96[tt],",
573 "[TT].pdb[tt] and [TT].g87[tt].",
574 "The file formats are detected from the file extension.",
575 "The precision of [TT].xtc[tt] and [TT].gro[tt] output is taken from the",
576 "input file for [TT].xtc[tt], [TT].gro[tt] and [TT].pdb[tt],",
577 "and from the [TT]-ndec[tt] option for other input formats. The precision",
578 "is always taken from [TT]-ndec[tt], when this option is set.",
579 "All other formats have fixed precision. [TT].trr[tt] and [TT].trj[tt]",
580 "output can be single or double precision, depending on the precision",
581 "of the [TT]trjconv[tt] binary.",
582 "Note that velocities are only supported in",
583 "[TT].trr[tt], [TT].trj[tt], [TT].gro[tt] and [TT].g96[tt] files.[PAR]",
585 "Option [TT]-app[tt] can be used to",
586 "append output to an existing trajectory file.",
587 "No checks are performed to ensure integrity",
588 "of the resulting combined trajectory file.[PAR]",
590 "Option [TT]-sep[tt] can be used to write every frame to a separate",
591 "[TT].gro, .g96[tt] or [TT].pdb[tt] file. By default, all frames all written to one file.",
592 "[TT].pdb[tt] files with all frames concatenated can be viewed with",
593 "[TT]rasmol -nmrpdb[tt].[PAR]",
595 "It is possible to select part of your trajectory and write it out",
596 "to a new trajectory file in order to save disk space, e.g. for leaving",
597 "out the water from a trajectory of a protein in water.",
598 "[BB]ALWAYS[bb] put the original trajectory on tape!",
599 "We recommend to use the portable [TT].xtc[tt] format for your analysis",
600 "to save disk space and to have portable files.[PAR]",
602 "There are two options for fitting the trajectory to a reference",
603 "either for essential dynamics analysis, etc.",
604 "The first option is just plain fitting to a reference structure",
605 "in the structure file. The second option is a progressive fit",
606 "in which the first timeframe is fitted to the reference structure ",
607 "in the structure file to obtain and each subsequent timeframe is ",
608 "fitted to the previously fitted structure. This way a continuous",
609 "trajectory is generated, which might not be the case when using the",
610 "regular fit method, e.g. when your protein undergoes large",
611 "conformational transitions.[PAR]",
613 "Option [TT]-pbc[tt] sets the type of periodic boundary condition",
615 "[TT]* mol[tt] puts the center of mass of molecules in the box,",
616 "and requires a run input file to be supplied with [TT]-s[tt].[BR]",
617 "[TT]* res[tt] puts the center of mass of residues in the box.[BR]",
618 "[TT]* atom[tt] puts all the atoms in the box.[BR]",
619 "[TT]* nojump[tt] checks if atoms jump across the box and then puts",
620 "them back. This has the effect that all molecules",
621 "will remain whole (provided they were whole in the initial",
622 "conformation). [BB]Note[bb] that this ensures a continuous trajectory but",
623 "molecules may diffuse out of the box. The starting configuration",
624 "for this procedure is taken from the structure file, if one is",
625 "supplied, otherwise it is the first frame.[BR]",
626 "[TT]* cluster[tt] clusters all the atoms in the selected index",
627 "such that they are all closest to the center of mass of the cluster,",
628 "which is iteratively updated. [BB]Note[bb] that this will only give meaningful",
629 "results if you in fact have a cluster. Luckily that can be checked",
630 "afterwards using a trajectory viewer. Note also that if your molecules",
631 "are broken this will not work either.[BR]",
632 "The separate option [TT]-clustercenter[tt] can be used to specify an",
633 "approximate center for the cluster. This is useful e.g. if you have",
634 "two big vesicles, and you want to maintain their relative positions.[BR]",
635 "[TT]* whole[tt] only makes broken molecules whole.[PAR]",
637 "Option [TT]-ur[tt] sets the unit cell representation for options",
638 "[TT]mol[tt], [TT]res[tt] and [TT]atom[tt] of [TT]-pbc[tt].",
639 "All three options give different results for triclinic boxes and",
640 "identical results for rectangular boxes.",
641 "[TT]rect[tt] is the ordinary brick shape.",
642 "[TT]tric[tt] is the triclinic unit cell.",
643 "[TT]compact[tt] puts all atoms at the closest distance from the center",
644 "of the box. This can be useful for visualizing e.g. truncated octahedra",
645 "or rhombic dodecahedra. The center for options [TT]tric[tt] and [TT]compact[tt]",
646 "is [TT]tric[tt] (see below), unless the option [TT]-boxcenter[tt]",
647 "is set differently.[PAR]",
649 "Option [TT]-center[tt] centers the system in the box. The user can",
650 "select the group which is used to determine the geometrical center.",
651 "Option [TT]-boxcenter[tt] sets the location of the center of the box",
652 "for options [TT]-pbc[tt] and [TT]-center[tt]. The center options are:",
653 "[TT]tric[tt]: half of the sum of the box vectors,",
654 "[TT]rect[tt]: half of the box diagonal,",
655 "[TT]zero[tt]: zero.",
656 "Use option [TT]-pbc mol[tt] in addition to [TT]-center[tt] when you",
657 "want all molecules in the box after the centering.[PAR]",
659 "It is not always possible to use combinations of [TT]-pbc[tt],",
660 "[TT]-fit[tt], [TT]-ur[tt] and [TT]-center[tt] to do exactly what",
661 "you want in one call to [TT]trjconv[tt]. Consider using multiple",
662 "calls, and check out the GROMACS website for suggestions.[PAR]",
664 "With [TT]-dt[tt], it is possible to reduce the number of ",
665 "frames in the output. This option relies on the accuracy of the times",
666 "in your input trajectory, so if these are inaccurate use the",
667 "[TT]-timestep[tt] option to modify the time (this can be done",
668 "simultaneously). For making smooth movies, the program [TT]g_filter[tt]",
669 "can reduce the number of frames while using low-pass frequency",
670 "filtering, this reduces aliasing of high frequency motions.[PAR]",
672 "Using [TT]-trunc[tt] [TT]trjconv[tt] can truncate [TT].trj[tt] in place, i.e.",
673 "without copying the file. This is useful when a run has crashed",
674 "during disk I/O (i.e. full disk), or when two contiguous",
675 "trajectories must be concatenated without having double frames.[PAR]",
677 "Option [TT]-dump[tt] can be used to extract a frame at or near",
678 "one specific time from your trajectory.[PAR]",
680 "Option [TT]-drop[tt] reads an [TT].xvg[tt] file with times and values.",
681 "When options [TT]-dropunder[tt] and/or [TT]-dropover[tt] are set,",
682 "frames with a value below and above the value of the respective options",
683 "will not be written."
699 const char *pbc_opt[epNR + 1] =
701 NULL, "none", "mol", "res", "atom", "nojump", "cluster", "whole",
706 const char *unitcell_opt[euNR+1] =
707 { NULL, "rect", "tric", "compact", NULL };
711 ecSel, ecTric, ecRect, ecZero, ecNR
713 const char *center_opt[ecNR+1] =
714 { NULL, "tric", "rect", "zero", NULL };
720 efSel, efNone, efFit, efFitXY, efReset, efResetXY, efPFit, efNR
722 const char *fit[efNR + 1] =
724 NULL, "none", "rot+trans", "rotxy+transxy", "translation", "transxy",
728 static gmx_bool bAppend = FALSE, bSeparate = FALSE, bVels = TRUE, bForce = FALSE, bCONECT = FALSE;
729 static gmx_bool bCenter = FALSE;
730 static int skip_nr = 1, ndec = 3, nzero = 0;
731 static real tzero = 0, delta_t = 0, timestep = 0, ttrunc = -1, tdump = -1, split_t = 0;
732 static rvec newbox = {0, 0, 0}, shift = {0, 0, 0}, trans = {0, 0, 0};
733 static char *exec_command = NULL;
734 static real dropunder = 0, dropover = 0;
735 static gmx_bool bRound = FALSE;
740 { "-skip", FALSE, etINT,
741 { &skip_nr }, "Only write every nr-th frame" },
742 { "-dt", FALSE, etTIME,
744 "Only write frame when t MOD dt = first time (%t)" },
745 { "-round", FALSE, etBOOL,
746 { &bRound }, "Round measurements to nearest picosecond"},
747 { "-dump", FALSE, etTIME,
748 { &tdump }, "Dump frame nearest specified time (%t)" },
749 { "-t0", FALSE, etTIME,
751 "Starting time (%t) (default: don't change)" },
752 { "-timestep", FALSE, etTIME,
754 "Change time step between input frames (%t)" },
755 { "-pbc", FALSE, etENUM,
757 "PBC treatment (see help text for full description)" },
758 { "-ur", FALSE, etENUM,
759 { unitcell_opt }, "Unit-cell representation" },
760 { "-center", FALSE, etBOOL,
761 { &bCenter }, "Center atoms in box" },
762 { "-boxcenter", FALSE, etENUM,
763 { center_opt }, "Center for -pbc and -center" },
764 { "-box", FALSE, etRVEC,
766 "Size for new cubic box (default: read from input)" },
767 { "-trans", FALSE, etRVEC,
769 "All coordinates will be translated by trans. This "
770 "can advantageously be combined with -pbc mol -ur "
772 { "-shift", FALSE, etRVEC,
774 "All coordinates will be shifted by framenr*shift" },
775 { "-fit", FALSE, etENUM,
777 "Fit molecule to ref structure in the structure file" },
778 { "-ndec", FALSE, etINT,
780 "Precision for .xtc and .gro writing in number of "
782 { "-vel", FALSE, etBOOL,
783 { &bVels }, "Read and write velocities if possible" },
784 { "-force", FALSE, etBOOL,
785 { &bForce }, "Read and write forces if possible" },
786 #ifndef GMX_NATIVE_WINDOWS
787 { "-trunc", FALSE, etTIME,
789 "Truncate input trajectory file after this time (%t)" },
791 { "-exec", FALSE, etSTR,
793 "Execute command for every output frame with the "
794 "frame number as argument" },
795 { "-app", FALSE, etBOOL,
796 { &bAppend }, "Append output" },
797 { "-split", FALSE, etTIME,
799 "Start writing new file when t MOD split = first "
801 { "-sep", FALSE, etBOOL,
803 "Write each frame to a separate .gro, .g96 or .pdb "
805 { "-nzero", FALSE, etINT,
807 "If the -sep flag is set, use these many digits "
808 "for the file numbers and prepend zeros as needed" },
809 { "-dropunder", FALSE, etREAL,
810 { &dropunder }, "Drop all frames below this value" },
811 { "-dropover", FALSE, etREAL,
812 { &dropover }, "Drop all frames above this value" },
813 { "-conect", FALSE, etBOOL,
815 "Add conect records when writing [TT].pdb[tt] files. Useful "
816 "for visualization of non-standard molecules, e.g. "
817 "coarse grained ones" }
819 #define NPA asize(pa)
822 t_trxstatus *trxout = NULL;
824 int ftp, ftpin = 0, file_nr;
825 t_trxframe fr, frout;
827 rvec *xmem = NULL, *vmem = NULL, *fmem = NULL;
828 rvec *xp = NULL, x_shift, hbox, box_center, dx;
829 real xtcpr, lambda, *w_rls = NULL;
830 int m, i, d, frame, outframe, natoms, nout, ncent, nre, newstep = 0, model_nr;
833 gmx_conect gc = NULL;
835 t_atoms *atoms = NULL, useatoms;
837 atom_id *index, *cindex;
841 int ifit, irms, my_clust = -1;
842 atom_id *ind_fit, *ind_rms;
843 char *gn_fit, *gn_rms;
844 t_cluster_ndx *clust = NULL;
845 t_trxstatus **clust_status = NULL;
846 int *clust_status_id = NULL;
848 int *nfwritten = NULL;
849 int ndrop = 0, ncol, drop0 = 0, drop1 = 0, dropuse = 0;
851 real tshift = 0, t0 = -1, dt = 0.001, prec;
852 gmx_bool bFit, bFitXY, bPFit, bReset;
854 gmx_rmpbc_t gpbc = NULL;
855 gmx_bool bRmPBC, bPBCWhole, bPBCcomRes, bPBCcomMol, bPBCcomAtom, bPBC, bNoJump, bCluster;
856 gmx_bool bCopy, bDoIt, bIndex, bTDump, bSetTime, bTPS = FALSE, bDTset = FALSE;
857 gmx_bool bExec, bTimeStep = FALSE, bDumpFrame = FALSE, bSetPrec, bNeedPrec;
858 gmx_bool bHaveFirstFrame, bHaveNextFrame, bSetBox, bSetUR, bSplit = FALSE;
859 gmx_bool bSubTraj = FALSE, bDropUnder = FALSE, bDropOver = FALSE, bTrans = FALSE;
860 gmx_bool bWriteFrame, bSplitHere;
861 const char *top_file, *in_file, *out_file = NULL;
862 char out_file2[256], *charpt;
863 char *outf_base = NULL;
864 const char *outf_ext = NULL;
865 char top_title[256], title[256], command[256], filemode[5];
867 gmx_bool bWarnCompact = FALSE;
872 { efTRX, "-f", NULL, ffREAD },
873 { efTRO, "-o", NULL, ffWRITE },
874 { efTPS, NULL, NULL, ffOPTRD },
875 { efNDX, NULL, NULL, ffOPTRD },
876 { efNDX, "-fr", "frames", ffOPTRD },
877 { efNDX, "-sub", "cluster", ffOPTRD },
878 { efXVG, "-drop", "drop", ffOPTRD }
880 #define NFILE asize(fnm)
882 parse_common_args(&argc, argv,
883 PCA_CAN_BEGIN | PCA_CAN_END | PCA_CAN_VIEW |
884 PCA_TIME_UNIT | PCA_BE_NICE,
885 NFILE, fnm, NPA, pa, asize(desc), desc,
888 top_file = ftp2fn(efTPS, NFILE, fnm);
891 /* Check command line */
892 in_file = opt2fn("-f", NFILE, fnm);
896 #ifndef GMX_NATIVE_WINDOWS
897 do_trunc(in_file, ttrunc);
902 /* mark active cmdline options */
903 bSetBox = opt2parg_bSet("-box", NPA, pa);
904 bSetTime = opt2parg_bSet("-t0", NPA, pa);
905 bSetPrec = opt2parg_bSet("-ndec", NPA, pa);
906 bSetUR = opt2parg_bSet("-ur", NPA, pa);
907 bExec = opt2parg_bSet("-exec", NPA, pa);
908 bTimeStep = opt2parg_bSet("-timestep", NPA, pa);
909 bTDump = opt2parg_bSet("-dump", NPA, pa);
910 bDropUnder = opt2parg_bSet("-dropunder", NPA, pa);
911 bDropOver = opt2parg_bSet("-dropover", NPA, pa);
912 bTrans = opt2parg_bSet("-trans", NPA, pa);
913 bSplit = (split_t != 0);
915 /* parse enum options */
916 fit_enum = nenum(fit);
917 bFit = (fit_enum == efFit || fit_enum == efFitXY);
918 bFitXY = fit_enum == efFitXY;
919 bReset = (fit_enum == efReset || fit_enum == efResetXY);
920 bPFit = fit_enum == efPFit;
921 pbc_enum = nenum(pbc_opt);
922 bPBCWhole = pbc_enum == epWhole;
923 bPBCcomRes = pbc_enum == epComRes;
924 bPBCcomMol = pbc_enum == epComMol;
925 bPBCcomAtom = pbc_enum == epComAtom;
926 bNoJump = pbc_enum == epNojump;
927 bCluster = pbc_enum == epCluster;
928 bPBC = pbc_enum != epNone;
929 unitcell_enum = nenum(unitcell_opt);
930 ecenter = nenum(center_opt) - ecTric;
932 /* set and check option dependencies */
935 bFit = TRUE; /* for pfit, fit *must* be set */
939 bReset = TRUE; /* for fit, reset *must* be set */
944 nfitdim = (fit_enum == efFitXY || fit_enum == efResetXY) ? 2 : 3;
946 bRmPBC = bFit || bPBCWhole || bPBCcomRes || bPBCcomMol;
950 if (!(bPBCcomRes || bPBCcomMol || bPBCcomAtom))
953 "WARNING: Option for unitcell representation (-ur %s)\n"
954 " only has effect in combination with -pbc %s, %s or %s.\n"
955 " Ingoring unitcell representation.\n\n",
956 unitcell_opt[0], pbc_opt[2], pbc_opt[3], pbc_opt[4]);
962 gmx_fatal(FARGS, "PBC condition treatment does not work together with rotational fit.\n"
963 "Please do the PBC condition treatment first and then run trjconv in a second step\n"
964 "for the rotational fit.\n"
965 "First doing the rotational fit and then doing the PBC treatment gives incorrect\n"
969 /* ndec is in nr of decimal places, prec is a multiplication factor: */
971 for (i = 0; i < ndec; i++)
976 bIndex = ftp2bSet(efNDX, NFILE, fnm);
979 /* Determine output type */
980 out_file = opt2fn("-o", NFILE, fnm);
981 ftp = fn2ftp(out_file);
982 fprintf(stderr, "Will write %s: %s\n", ftp2ext(ftp), ftp2desc(ftp));
983 bNeedPrec = (ftp == efXTC || ftp == efGRO);
986 /* check if velocities are possible in input and output files */
987 ftpin = fn2ftp(in_file);
988 bVels = (ftp == efTRR || ftp == efTRJ || ftp == efGRO || ftp == efG96)
989 && (ftpin == efTRR || ftpin == efTRJ || ftpin == efGRO || ftpin == efG96 ||
992 if (bSeparate || bSplit)
994 outf_ext = strrchr(out_file, '.');
995 if (outf_ext == NULL)
997 gmx_fatal(FARGS, "Output file name '%s' does not contain a '.'", out_file);
999 outf_base = strdup(out_file);
1000 outf_base[outf_ext - out_file] = '\0';
1003 bSubTraj = opt2bSet("-sub", NFILE, fnm);
1006 if ((ftp != efXTC) && (ftp != efTRR))
1008 /* It seems likely that other trajectory file types
1009 * could work here. */
1010 gmx_fatal(FARGS, "Can only use the sub option with output file types "
1013 clust = cluster_index(NULL, opt2fn("-sub", NFILE, fnm));
1015 /* Check for number of files disabled, as FOPEN_MAX is not the correct
1016 * number to check for. In my linux box it is only 16.
1018 if (0 && (clust->clust->nr > FOPEN_MAX-4))
1020 gmx_fatal(FARGS, "Can not open enough (%d) files to write all the"
1021 " trajectories.\ntry splitting the index file in %d parts.\n"
1023 clust->clust->nr, 1+clust->clust->nr/FOPEN_MAX, FOPEN_MAX);
1025 gmx_warning("The -sub option could require as many open output files as there are\n"
1026 "index groups in the file (%d). If you get I/O errors opening new files,\n"
1027 "try reducing the number of index groups in the file, and perhaps\n"
1028 "using trjconv -sub several times on different chunks of your index file.\n",
1031 snew(clust_status, clust->clust->nr);
1032 snew(clust_status_id, clust->clust->nr);
1033 snew(nfwritten, clust->clust->nr);
1034 for (i = 0; (i < clust->clust->nr); i++)
1036 clust_status[i] = NULL;
1037 clust_status_id[i] = -1;
1039 bSeparate = bSplit = FALSE;
1046 /* Determine whether to read a topology */
1047 bTPS = (ftp2bSet(efTPS, NFILE, fnm) ||
1048 bRmPBC || bReset || bPBCcomMol || bCluster ||
1049 (ftp == efGRO) || (ftp == efPDB) || bCONECT);
1051 /* Determine if when can read index groups */
1052 bIndex = (bIndex || bTPS);
1056 read_tps_conf(top_file, top_title, &top, &ePBC, &xp, NULL, top_box,
1057 bReset || bPBCcomRes);
1060 if (0 == top.mols.nr && (bCluster || bPBCcomMol))
1062 gmx_fatal(FARGS, "Option -pbc %s requires a .tpr file for the -s option", pbc_opt[pbc_enum]);
1065 /* top_title is only used for gro and pdb,
1066 * the header in such a file is top_title t= ...
1067 * to prevent a double t=, remove it from top_title
1069 if ((charpt = strstr(top_title, " t= ")))
1076 gc = gmx_conect_generate(&top);
1080 gpbc = gmx_rmpbc_init(&top.idef, ePBC, top.atoms.nr);
1084 /* get frame number index */
1086 if (opt2bSet("-fr", NFILE, fnm))
1088 printf("Select groups of frame number indices:\n");
1089 rd_index(opt2fn("-fr", NFILE, fnm), 1, &nrfri, (atom_id **)&frindex, &frname);
1092 for (i = 0; i < nrfri; i++)
1094 fprintf(debug, "frindex[%4d]=%4d\n", i, frindex[i]);
1099 /* get index groups etc. */
1102 printf("Select group for %s fit\n",
1103 bFit ? "least squares" : "translational");
1104 get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm),
1105 1, &ifit, &ind_fit, &gn_fit);
1111 gmx_fatal(FARGS, "Need at least 2 atoms to fit!\n");
1115 fprintf(stderr, "WARNING: fitting with only 2 atoms is not unique\n");
1121 printf("Select group for clustering\n");
1122 get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm),
1123 1, &ifit, &ind_fit, &gn_fit);
1130 printf("Select group for centering\n");
1131 get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm),
1132 1, &ncent, &cindex, &grpnm);
1134 printf("Select group for output\n");
1135 get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm),
1136 1, &nout, &index, &grpnm);
1140 /* no index file, so read natoms from TRX */
1141 if (!read_first_frame(oenv, &status, in_file, &fr, TRX_DONT_SKIP))
1143 gmx_fatal(FARGS, "Could not read a frame from %s", in_file);
1148 snew(index, natoms);
1149 for (i = 0; i < natoms; i++)
1163 snew(w_rls, atoms->nr);
1164 for (i = 0; (i < ifit); i++)
1166 w_rls[ind_fit[i]] = atoms->atom[ind_fit[i]].m;
1169 /* Restore reference structure and set to origin,
1170 store original location (to put structure back) */
1173 gmx_rmpbc(gpbc, top.atoms.nr, top_box, xp);
1175 copy_rvec(xp[index[0]], x_shift);
1176 reset_x_ndim(nfitdim, ifit, ind_fit, atoms->nr, NULL, xp, w_rls);
1177 rvec_dec(x_shift, xp[index[0]]);
1181 clear_rvec(x_shift);
1184 if (bDropUnder || bDropOver)
1186 /* Read the .xvg file with the drop values */
1187 fprintf(stderr, "\nReading drop file ...");
1188 ndrop = read_xvg(opt2fn("-drop", NFILE, fnm), &dropval, &ncol);
1189 fprintf(stderr, " %d time points\n", ndrop);
1190 if (ndrop == 0 || ncol < 2)
1192 gmx_fatal(FARGS, "Found no data points in %s",
1193 opt2fn("-drop", NFILE, fnm));
1199 /* Make atoms struct for output in GRO or PDB files */
1200 if ((ftp == efGRO) || ((ftp == efG96) && bTPS) || (ftp == efPDB))
1202 /* get memory for stuff to go in .pdb file */
1203 init_t_atoms(&useatoms, atoms->nr, FALSE);
1204 sfree(useatoms.resinfo);
1205 useatoms.resinfo = atoms->resinfo;
1206 for (i = 0; (i < nout); i++)
1208 useatoms.atomname[i] = atoms->atomname[index[i]];
1209 useatoms.atom[i] = atoms->atom[index[i]];
1210 useatoms.nres = max(useatoms.nres, useatoms.atom[i].resind+1);
1214 /* select what to read */
1215 if (ftp == efTRR || ftp == efTRJ)
1225 flags = flags | TRX_READ_V;
1229 flags = flags | TRX_READ_F;
1232 /* open trx file for reading */
1233 bHaveFirstFrame = read_first_frame(oenv, &status, in_file, &fr, flags);
1236 fprintf(stderr, "\nPrecision of %s is %g (nm)\n", in_file, 1/fr.prec);
1240 if (bSetPrec || !fr.bPrec)
1242 fprintf(stderr, "\nSetting output precision to %g (nm)\n", 1/prec);
1246 fprintf(stderr, "Using output precision of %g (nm)\n", 1/prec);
1250 if (bHaveFirstFrame)
1252 set_trxframe_ePBC(&fr, ePBC);
1258 tshift = tzero-fr.time;
1265 /* open output for writing */
1266 if ((bAppend) && (gmx_fexist(out_file)))
1268 strcpy(filemode, "a");
1269 fprintf(stderr, "APPENDING to existing file %s\n", out_file);
1273 strcpy(filemode, "w");
1282 if (!bSplit && !bSubTraj)
1284 trxout = open_trx(out_file, filemode);
1290 if (( !bSeparate && !bSplit ) && !bSubTraj)
1292 out = ffopen(out_file, filemode);
1300 /* check if index is meaningful */
1301 for (i = 0; i < nout; i++)
1303 if (index[i] >= natoms)
1306 "Index[%d] %d is larger than the number of atoms in the\n"
1307 "trajectory file (%d). There is a mismatch in the contents\n"
1308 "of your -f, -s and/or -n files.", i, index[i]+1, natoms);
1310 bCopy = bCopy || (i != index[i]);
1328 fprintf(gmx_fio_getfp(trx_get_fileio(trxout)),
1329 "Generated by %s. #atoms=%d, a BOX is"
1330 " stored in this file.\n", Program(), nout);
1333 /* Start the big loop over frames */
1340 /* Main loop over frames */
1351 /*if (frame >= clust->clust->nra)
1352 gmx_fatal(FARGS,"There are more frames in the trajectory than in the cluster index file\n");*/
1353 if (frame > clust->maxframe)
1359 my_clust = clust->inv_clust[frame];
1361 if ((my_clust < 0) || (my_clust >= clust->clust->nr) ||
1362 (my_clust == NO_ATID))
1370 /* generate new box */
1372 for (m = 0; m < DIM; m++)
1374 fr.box[m][m] = newbox[m];
1380 for (i = 0; i < natoms; i++)
1382 rvec_inc(fr.x[i], trans);
1388 /* determine timestep */
1401 /* This is not very elegant, as one can not dump a frame after
1402 * a timestep with is more than twice as small as the first one. */
1403 bDumpFrame = (fr.time > tdump-0.5*dt) && (fr.time <= tdump+0.5*dt);
1410 /* determine if an atom jumped across the box and reset it if so */
1411 if (bNoJump && (bTPS || frame != 0))
1413 for (d = 0; d < DIM; d++)
1415 hbox[d] = 0.5*fr.box[d][d];
1417 for (i = 0; i < natoms; i++)
1421 rvec_dec(fr.x[i], x_shift);
1423 for (m = DIM-1; m >= 0; m--)
1427 while (fr.x[i][m]-xp[i][m] <= -hbox[m])
1429 for (d = 0; d <= m; d++)
1431 fr.x[i][d] += fr.box[m][d];
1434 while (fr.x[i][m]-xp[i][m] > hbox[m])
1436 for (d = 0; d <= m; d++)
1438 fr.x[i][d] -= fr.box[m][d];
1447 calc_pbc_cluster(ecenter, ifit, &top, ePBC, fr.x, ind_fit, fr.box);
1452 /* Now modify the coords according to the flags,
1453 for normal fit, this is only done for output frames */
1456 gmx_rmpbc_trxfr(gpbc, &fr);
1459 reset_x_ndim(nfitdim, ifit, ind_fit, natoms, NULL, fr.x, w_rls);
1460 do_fit(natoms, w_rls, xp, fr.x);
1463 /* store this set of coordinates for future use */
1464 if (bPFit || bNoJump)
1470 for (i = 0; (i < natoms); i++)
1472 copy_rvec(fr.x[i], xp[i]);
1473 rvec_inc(fr.x[i], x_shift);
1479 /* see if we have a frame from the frame index group */
1480 for (i = 0; i < nrfri && !bDumpFrame; i++)
1482 bDumpFrame = frame == frindex[i];
1485 if (debug && bDumpFrame)
1487 fprintf(debug, "dumping %d\n", frame);
1491 ( ( !bTDump && !frindex && frame % skip_nr == 0 ) || bDumpFrame );
1493 if (bWriteFrame && (bDropUnder || bDropOver))
1495 while (dropval[0][drop1] < fr.time && drop1+1 < ndrop)
1500 if (fabs(dropval[0][drop0] - fr.time)
1501 < fabs(dropval[0][drop1] - fr.time))
1509 if ((bDropUnder && dropval[1][dropuse] < dropunder) ||
1510 (bDropOver && dropval[1][dropuse] > dropover))
1512 bWriteFrame = FALSE;
1522 fr.time = tzero+frame*timestep;
1532 fprintf(stderr, "\nDumping frame at t= %g %s\n",
1533 output_env_conv_time(oenv, fr.time), output_env_get_time_unit(oenv));
1536 /* check for writing at each delta_t */
1537 bDoIt = (delta_t == 0);
1542 bDoIt = bRmod(fr.time, tzero, delta_t);
1546 /* round() is not C89 compatible, so we do this: */
1547 bDoIt = bRmod(floor(fr.time+0.5), floor(tzero+0.5),
1548 floor(delta_t+0.5));
1552 if (bDoIt || bTDump)
1554 /* print sometimes */
1555 if ( ((outframe % SKIP) == 0) || (outframe < SKIP) )
1557 fprintf(stderr, " -> frame %6d time %8.3f \r",
1558 outframe, output_env_conv_time(oenv, fr.time));
1563 /* Now modify the coords according to the flags,
1564 for PFit we did this already! */
1568 gmx_rmpbc_trxfr(gpbc, &fr);
1573 reset_x_ndim(nfitdim, ifit, ind_fit, natoms, NULL, fr.x, w_rls);
1576 do_fit_ndim(nfitdim, natoms, w_rls, xp, fr.x);
1580 for (i = 0; i < natoms; i++)
1582 rvec_inc(fr.x[i], x_shift);
1589 center_x(ecenter, fr.x, fr.box, natoms, ncent, cindex);
1595 switch (unitcell_enum)
1598 put_atoms_in_box(ePBC, fr.box, natoms, fr.x);
1601 put_atoms_in_triclinic_unitcell(ecenter, fr.box, natoms, fr.x);
1604 warn = put_atoms_in_compact_unitcell(ePBC, ecenter, fr.box,
1606 if (warn && !bWarnCompact)
1608 fprintf(stderr, "\n%s\n", warn);
1609 bWarnCompact = TRUE;
1616 put_residue_com_in_box(unitcell_enum, ecenter,
1617 natoms, atoms->atom, ePBC, fr.box, fr.x);
1621 put_molecule_com_in_box(unitcell_enum, ecenter,
1623 natoms, atoms->atom, ePBC, fr.box, fr.x);
1625 /* Copy the input trxframe struct to the output trxframe struct */
1627 frout.bV = (frout.bV && bVels);
1628 frout.bF = (frout.bF && bForce);
1629 frout.natoms = nout;
1630 if (bNeedPrec && (bSetPrec || !fr.bPrec))
1646 for (i = 0; i < nout; i++)
1648 copy_rvec(fr.x[index[i]], frout.x[i]);
1651 copy_rvec(fr.v[index[i]], frout.v[i]);
1655 copy_rvec(fr.f[index[i]], frout.f[i]);
1660 if (opt2parg_bSet("-shift", NPA, pa))
1662 for (i = 0; i < nout; i++)
1664 for (d = 0; d < DIM; d++)
1666 frout.x[i][d] += outframe*shift[d];
1673 bSplitHere = bSplit && bRmod(fr.time, tzero, split_t);
1677 /* round() is not C89 compatible, so we do this: */
1678 bSplitHere = bSplit && bRmod(floor(fr.time+0.5),
1680 floor(split_t+0.5));
1682 if (bSeparate || bSplitHere)
1684 mk_filenm(outf_base, ftp2ext(ftp), nzero, file_nr, out_file2);
1699 trxout = open_trx(out_file2, filemode);
1706 if (clust_status_id[my_clust] == -1)
1708 sprintf(buf, "%s.%s", clust->grpname[my_clust], ftp2ext(ftp));
1709 clust_status[my_clust] = open_trx(buf, "w");
1710 clust_status_id[my_clust] = 1;
1713 else if (clust_status_id[my_clust] == -2)
1715 gmx_fatal(FARGS, "File %s.xtc should still be open (%d open .xtc files)\n" "in order to write frame %d. my_clust = %d",
1716 clust->grpname[my_clust], ntrxopen, frame,
1719 write_trxframe(clust_status[my_clust], &frout, gc);
1720 nfwritten[my_clust]++;
1721 if (nfwritten[my_clust] ==
1722 (clust->clust->index[my_clust+1]-
1723 clust->clust->index[my_clust]))
1725 close_trx(clust_status[my_clust]);
1726 clust_status[my_clust] = NULL;
1727 clust_status_id[my_clust] = -2;
1731 gmx_fatal(FARGS, "Less than zero open .xtc files!");
1738 write_trxframe(trxout, &frout, gc);
1744 sprintf(title, "Generated by trjconv : %s t= %9.5f",
1745 top_title, fr.time);
1746 if (bSeparate || bSplitHere)
1748 out = ffopen(out_file2, "w");
1753 write_hconf_p(out, title, &useatoms, prec2ndec(frout.prec),
1754 frout.x, frout.bV ? frout.v : NULL, frout.box);
1757 fprintf(out, "REMARK GENERATED BY TRJCONV\n");
1758 sprintf(title, "%s t= %9.5f", top_title, frout.time);
1759 /* if reading from pdb, we want to keep the original
1760 model numbering else we write the output frame
1761 number plus one, because model 0 is not allowed in pdb */
1762 if (ftpin == efPDB && fr.bStep && fr.step > model_nr)
1770 write_pdbfile(out, title, &useatoms, frout.x,
1771 frout.ePBC, frout.box, ' ', model_nr, gc, TRUE);
1774 frout.title = title;
1775 if (bSeparate || bTDump)
1777 frout.bTitle = TRUE;
1780 frout.bAtoms = TRUE;
1782 frout.atoms = &useatoms;
1783 frout.bStep = FALSE;
1784 frout.bTime = FALSE;
1788 frout.bTitle = (outframe == 0);
1789 frout.bAtoms = FALSE;
1793 write_g96_conf(out, &frout, -1, NULL);
1802 gmx_fatal(FARGS, "DHE, ftp=%d\n", ftp);
1804 if (bSeparate || bSplitHere)
1809 /* execute command */
1813 sprintf(c, "%s %d", exec_command, file_nr-1);
1814 /*fprintf(stderr,"Executing '%s'\n",c);*/
1815 #ifdef GMX_NO_SYSTEM
1816 printf("Warning-- No calls to system(3) supported on this platform.");
1817 printf("Warning-- Skipping execution of 'system(\"%s\")'.", c);
1821 gmx_fatal(FARGS, "Error executing command: %s", c);
1829 bHaveNextFrame = read_next_frame(oenv, status, &fr);
1831 while (!(bTDump && bDumpFrame) && bHaveNextFrame);
1834 if (!bHaveFirstFrame || (bTDump && !bDumpFrame))
1836 fprintf(stderr, "\nWARNING no output, "
1837 "last frame read at t=%g\n", fr.time);
1839 fprintf(stderr, "\n");
1846 gmx_rmpbc_done(gpbc);
1853 else if (out != NULL)
1859 for (i = 0; (i < clust->clust->nr); i++)
1861 if (clust_status_id[i] >= 0)
1863 close_trx(clust_status[i]);
1869 do_view(oenv, out_file, NULL);