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 if (!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,
891 top_file = ftp2fn(efTPS, NFILE, fnm);
894 /* Check command line */
895 in_file = opt2fn("-f", NFILE, fnm);
899 #ifndef GMX_NATIVE_WINDOWS
900 do_trunc(in_file, ttrunc);
905 /* mark active cmdline options */
906 bSetBox = opt2parg_bSet("-box", NPA, pa);
907 bSetTime = opt2parg_bSet("-t0", NPA, pa);
908 bSetPrec = opt2parg_bSet("-ndec", NPA, pa);
909 bSetUR = opt2parg_bSet("-ur", NPA, pa);
910 bExec = opt2parg_bSet("-exec", NPA, pa);
911 bTimeStep = opt2parg_bSet("-timestep", NPA, pa);
912 bTDump = opt2parg_bSet("-dump", NPA, pa);
913 bDropUnder = opt2parg_bSet("-dropunder", NPA, pa);
914 bDropOver = opt2parg_bSet("-dropover", NPA, pa);
915 bTrans = opt2parg_bSet("-trans", NPA, pa);
916 bSplit = (split_t != 0);
918 /* parse enum options */
919 fit_enum = nenum(fit);
920 bFit = (fit_enum == efFit || fit_enum == efFitXY);
921 bFitXY = fit_enum == efFitXY;
922 bReset = (fit_enum == efReset || fit_enum == efResetXY);
923 bPFit = fit_enum == efPFit;
924 pbc_enum = nenum(pbc_opt);
925 bPBCWhole = pbc_enum == epWhole;
926 bPBCcomRes = pbc_enum == epComRes;
927 bPBCcomMol = pbc_enum == epComMol;
928 bPBCcomAtom = pbc_enum == epComAtom;
929 bNoJump = pbc_enum == epNojump;
930 bCluster = pbc_enum == epCluster;
931 bPBC = pbc_enum != epNone;
932 unitcell_enum = nenum(unitcell_opt);
933 ecenter = nenum(center_opt) - ecTric;
935 /* set and check option dependencies */
938 bFit = TRUE; /* for pfit, fit *must* be set */
942 bReset = TRUE; /* for fit, reset *must* be set */
947 nfitdim = (fit_enum == efFitXY || fit_enum == efResetXY) ? 2 : 3;
949 bRmPBC = bFit || bPBCWhole || bPBCcomRes || bPBCcomMol;
953 if (!(bPBCcomRes || bPBCcomMol || bPBCcomAtom))
956 "WARNING: Option for unitcell representation (-ur %s)\n"
957 " only has effect in combination with -pbc %s, %s or %s.\n"
958 " Ingoring unitcell representation.\n\n",
959 unitcell_opt[0], pbc_opt[2], pbc_opt[3], pbc_opt[4]);
965 gmx_fatal(FARGS, "PBC condition treatment does not work together with rotational fit.\n"
966 "Please do the PBC condition treatment first and then run trjconv in a second step\n"
967 "for the rotational fit.\n"
968 "First doing the rotational fit and then doing the PBC treatment gives incorrect\n"
972 /* ndec is in nr of decimal places, prec is a multiplication factor: */
974 for (i = 0; i < ndec; i++)
979 bIndex = ftp2bSet(efNDX, NFILE, fnm);
982 /* Determine output type */
983 out_file = opt2fn("-o", NFILE, fnm);
984 ftp = fn2ftp(out_file);
985 fprintf(stderr, "Will write %s: %s\n", ftp2ext(ftp), ftp2desc(ftp));
986 bNeedPrec = (ftp == efXTC || ftp == efGRO);
989 /* check if velocities are possible in input and output files */
990 ftpin = fn2ftp(in_file);
991 bVels = (ftp == efTRR || ftp == efTRJ || ftp == efGRO || ftp == efG96)
992 && (ftpin == efTRR || ftpin == efTRJ || ftpin == efGRO || ftpin == efG96 ||
995 if (bSeparate || bSplit)
997 outf_ext = strrchr(out_file, '.');
998 if (outf_ext == NULL)
1000 gmx_fatal(FARGS, "Output file name '%s' does not contain a '.'", out_file);
1002 outf_base = strdup(out_file);
1003 outf_base[outf_ext - out_file] = '\0';
1006 bSubTraj = opt2bSet("-sub", NFILE, fnm);
1009 if ((ftp != efXTC) && (ftp != efTRR))
1011 /* It seems likely that other trajectory file types
1012 * could work here. */
1013 gmx_fatal(FARGS, "Can only use the sub option with output file types "
1016 clust = cluster_index(NULL, opt2fn("-sub", NFILE, fnm));
1018 /* Check for number of files disabled, as FOPEN_MAX is not the correct
1019 * number to check for. In my linux box it is only 16.
1021 if (0 && (clust->clust->nr > FOPEN_MAX-4))
1023 gmx_fatal(FARGS, "Can not open enough (%d) files to write all the"
1024 " trajectories.\ntry splitting the index file in %d parts.\n"
1026 clust->clust->nr, 1+clust->clust->nr/FOPEN_MAX, FOPEN_MAX);
1028 gmx_warning("The -sub option could require as many open output files as there are\n"
1029 "index groups in the file (%d). If you get I/O errors opening new files,\n"
1030 "try reducing the number of index groups in the file, and perhaps\n"
1031 "using trjconv -sub several times on different chunks of your index file.\n",
1034 snew(clust_status, clust->clust->nr);
1035 snew(clust_status_id, clust->clust->nr);
1036 snew(nfwritten, clust->clust->nr);
1037 for (i = 0; (i < clust->clust->nr); i++)
1039 clust_status[i] = NULL;
1040 clust_status_id[i] = -1;
1042 bSeparate = bSplit = FALSE;
1049 /* Determine whether to read a topology */
1050 bTPS = (ftp2bSet(efTPS, NFILE, fnm) ||
1051 bRmPBC || bReset || bPBCcomMol || bCluster ||
1052 (ftp == efGRO) || (ftp == efPDB) || bCONECT);
1054 /* Determine if when can read index groups */
1055 bIndex = (bIndex || bTPS);
1059 read_tps_conf(top_file, top_title, &top, &ePBC, &xp, NULL, top_box,
1060 bReset || bPBCcomRes);
1063 if (0 == top.mols.nr && (bCluster || bPBCcomMol))
1065 gmx_fatal(FARGS, "Option -pbc %s requires a .tpr file for the -s option", pbc_opt[pbc_enum]);
1068 /* top_title is only used for gro and pdb,
1069 * the header in such a file is top_title t= ...
1070 * to prevent a double t=, remove it from top_title
1072 if ((charpt = strstr(top_title, " t= ")))
1079 gc = gmx_conect_generate(&top);
1083 gpbc = gmx_rmpbc_init(&top.idef, ePBC, top.atoms.nr);
1087 /* get frame number index */
1089 if (opt2bSet("-fr", NFILE, fnm))
1091 printf("Select groups of frame number indices:\n");
1092 rd_index(opt2fn("-fr", NFILE, fnm), 1, &nrfri, (atom_id **)&frindex, &frname);
1095 for (i = 0; i < nrfri; i++)
1097 fprintf(debug, "frindex[%4d]=%4d\n", i, frindex[i]);
1102 /* get index groups etc. */
1105 printf("Select group for %s fit\n",
1106 bFit ? "least squares" : "translational");
1107 get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm),
1108 1, &ifit, &ind_fit, &gn_fit);
1114 gmx_fatal(FARGS, "Need at least 2 atoms to fit!\n");
1118 fprintf(stderr, "WARNING: fitting with only 2 atoms is not unique\n");
1124 printf("Select group for clustering\n");
1125 get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm),
1126 1, &ifit, &ind_fit, &gn_fit);
1133 printf("Select group for centering\n");
1134 get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm),
1135 1, &ncent, &cindex, &grpnm);
1137 printf("Select group for output\n");
1138 get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm),
1139 1, &nout, &index, &grpnm);
1143 /* no index file, so read natoms from TRX */
1144 if (!read_first_frame(oenv, &status, in_file, &fr, TRX_DONT_SKIP))
1146 gmx_fatal(FARGS, "Could not read a frame from %s", in_file);
1151 snew(index, natoms);
1152 for (i = 0; i < natoms; i++)
1166 snew(w_rls, atoms->nr);
1167 for (i = 0; (i < ifit); i++)
1169 w_rls[ind_fit[i]] = atoms->atom[ind_fit[i]].m;
1172 /* Restore reference structure and set to origin,
1173 store original location (to put structure back) */
1176 gmx_rmpbc(gpbc, top.atoms.nr, top_box, xp);
1178 copy_rvec(xp[index[0]], x_shift);
1179 reset_x_ndim(nfitdim, ifit, ind_fit, atoms->nr, NULL, xp, w_rls);
1180 rvec_dec(x_shift, xp[index[0]]);
1184 clear_rvec(x_shift);
1187 if (bDropUnder || bDropOver)
1189 /* Read the .xvg file with the drop values */
1190 fprintf(stderr, "\nReading drop file ...");
1191 ndrop = read_xvg(opt2fn("-drop", NFILE, fnm), &dropval, &ncol);
1192 fprintf(stderr, " %d time points\n", ndrop);
1193 if (ndrop == 0 || ncol < 2)
1195 gmx_fatal(FARGS, "Found no data points in %s",
1196 opt2fn("-drop", NFILE, fnm));
1202 /* Make atoms struct for output in GRO or PDB files */
1203 if ((ftp == efGRO) || ((ftp == efG96) && bTPS) || (ftp == efPDB))
1205 /* get memory for stuff to go in .pdb file */
1206 init_t_atoms(&useatoms, atoms->nr, FALSE);
1207 sfree(useatoms.resinfo);
1208 useatoms.resinfo = atoms->resinfo;
1209 for (i = 0; (i < nout); i++)
1211 useatoms.atomname[i] = atoms->atomname[index[i]];
1212 useatoms.atom[i] = atoms->atom[index[i]];
1213 useatoms.nres = max(useatoms.nres, useatoms.atom[i].resind+1);
1217 /* select what to read */
1218 if (ftp == efTRR || ftp == efTRJ)
1228 flags = flags | TRX_READ_V;
1232 flags = flags | TRX_READ_F;
1235 /* open trx file for reading */
1236 bHaveFirstFrame = read_first_frame(oenv, &status, in_file, &fr, flags);
1239 fprintf(stderr, "\nPrecision of %s is %g (nm)\n", in_file, 1/fr.prec);
1243 if (bSetPrec || !fr.bPrec)
1245 fprintf(stderr, "\nSetting output precision to %g (nm)\n", 1/prec);
1249 fprintf(stderr, "Using output precision of %g (nm)\n", 1/prec);
1253 if (bHaveFirstFrame)
1255 set_trxframe_ePBC(&fr, ePBC);
1261 tshift = tzero-fr.time;
1268 /* open output for writing */
1269 if ((bAppend) && (gmx_fexist(out_file)))
1271 strcpy(filemode, "a");
1272 fprintf(stderr, "APPENDING to existing file %s\n", out_file);
1276 strcpy(filemode, "w");
1285 if (!bSplit && !bSubTraj)
1287 trxout = open_trx(out_file, filemode);
1293 if (( !bSeparate && !bSplit ) && !bSubTraj)
1295 out = ffopen(out_file, filemode);
1303 /* check if index is meaningful */
1304 for (i = 0; i < nout; i++)
1306 if (index[i] >= natoms)
1309 "Index[%d] %d is larger than the number of atoms in the\n"
1310 "trajectory file (%d). There is a mismatch in the contents\n"
1311 "of your -f, -s and/or -n files.", i, index[i]+1, natoms);
1313 bCopy = bCopy || (i != index[i]);
1331 fprintf(gmx_fio_getfp(trx_get_fileio(trxout)),
1332 "Generated by %s. #atoms=%d, a BOX is"
1333 " stored in this file.\n", ShortProgram(), nout);
1336 /* Start the big loop over frames */
1343 /* Main loop over frames */
1354 /*if (frame >= clust->clust->nra)
1355 gmx_fatal(FARGS,"There are more frames in the trajectory than in the cluster index file\n");*/
1356 if (frame > clust->maxframe)
1362 my_clust = clust->inv_clust[frame];
1364 if ((my_clust < 0) || (my_clust >= clust->clust->nr) ||
1365 (my_clust == NO_ATID))
1373 /* generate new box */
1375 for (m = 0; m < DIM; m++)
1377 fr.box[m][m] = newbox[m];
1383 for (i = 0; i < natoms; i++)
1385 rvec_inc(fr.x[i], trans);
1391 /* determine timestep */
1404 /* This is not very elegant, as one can not dump a frame after
1405 * a timestep with is more than twice as small as the first one. */
1406 bDumpFrame = (fr.time > tdump-0.5*dt) && (fr.time <= tdump+0.5*dt);
1413 /* determine if an atom jumped across the box and reset it if so */
1414 if (bNoJump && (bTPS || frame != 0))
1416 for (d = 0; d < DIM; d++)
1418 hbox[d] = 0.5*fr.box[d][d];
1420 for (i = 0; i < natoms; i++)
1424 rvec_dec(fr.x[i], x_shift);
1426 for (m = DIM-1; m >= 0; m--)
1430 while (fr.x[i][m]-xp[i][m] <= -hbox[m])
1432 for (d = 0; d <= m; d++)
1434 fr.x[i][d] += fr.box[m][d];
1437 while (fr.x[i][m]-xp[i][m] > hbox[m])
1439 for (d = 0; d <= m; d++)
1441 fr.x[i][d] -= fr.box[m][d];
1450 calc_pbc_cluster(ecenter, ifit, &top, ePBC, fr.x, ind_fit, fr.box);
1455 /* Now modify the coords according to the flags,
1456 for normal fit, this is only done for output frames */
1459 gmx_rmpbc_trxfr(gpbc, &fr);
1462 reset_x_ndim(nfitdim, ifit, ind_fit, natoms, NULL, fr.x, w_rls);
1463 do_fit(natoms, w_rls, xp, fr.x);
1466 /* store this set of coordinates for future use */
1467 if (bPFit || bNoJump)
1473 for (i = 0; (i < natoms); i++)
1475 copy_rvec(fr.x[i], xp[i]);
1476 rvec_inc(fr.x[i], x_shift);
1482 /* see if we have a frame from the frame index group */
1483 for (i = 0; i < nrfri && !bDumpFrame; i++)
1485 bDumpFrame = frame == frindex[i];
1488 if (debug && bDumpFrame)
1490 fprintf(debug, "dumping %d\n", frame);
1494 ( ( !bTDump && !frindex && frame % skip_nr == 0 ) || bDumpFrame );
1496 if (bWriteFrame && (bDropUnder || bDropOver))
1498 while (dropval[0][drop1] < fr.time && drop1+1 < ndrop)
1503 if (fabs(dropval[0][drop0] - fr.time)
1504 < fabs(dropval[0][drop1] - fr.time))
1512 if ((bDropUnder && dropval[1][dropuse] < dropunder) ||
1513 (bDropOver && dropval[1][dropuse] > dropover))
1515 bWriteFrame = FALSE;
1525 fr.time = tzero+frame*timestep;
1535 fprintf(stderr, "\nDumping frame at t= %g %s\n",
1536 output_env_conv_time(oenv, fr.time), output_env_get_time_unit(oenv));
1539 /* check for writing at each delta_t */
1540 bDoIt = (delta_t == 0);
1545 bDoIt = bRmod(fr.time, tzero, delta_t);
1549 /* round() is not C89 compatible, so we do this: */
1550 bDoIt = bRmod(floor(fr.time+0.5), floor(tzero+0.5),
1551 floor(delta_t+0.5));
1555 if (bDoIt || bTDump)
1557 /* print sometimes */
1558 if ( ((outframe % SKIP) == 0) || (outframe < SKIP) )
1560 fprintf(stderr, " -> frame %6d time %8.3f \r",
1561 outframe, output_env_conv_time(oenv, fr.time));
1566 /* Now modify the coords according to the flags,
1567 for PFit we did this already! */
1571 gmx_rmpbc_trxfr(gpbc, &fr);
1576 reset_x_ndim(nfitdim, ifit, ind_fit, natoms, NULL, fr.x, w_rls);
1579 do_fit_ndim(nfitdim, natoms, w_rls, xp, fr.x);
1583 for (i = 0; i < natoms; i++)
1585 rvec_inc(fr.x[i], x_shift);
1592 center_x(ecenter, fr.x, fr.box, natoms, ncent, cindex);
1598 switch (unitcell_enum)
1601 put_atoms_in_box(ePBC, fr.box, natoms, fr.x);
1604 put_atoms_in_triclinic_unitcell(ecenter, fr.box, natoms, fr.x);
1607 warn = put_atoms_in_compact_unitcell(ePBC, ecenter, fr.box,
1609 if (warn && !bWarnCompact)
1611 fprintf(stderr, "\n%s\n", warn);
1612 bWarnCompact = TRUE;
1619 put_residue_com_in_box(unitcell_enum, ecenter,
1620 natoms, atoms->atom, ePBC, fr.box, fr.x);
1624 put_molecule_com_in_box(unitcell_enum, ecenter,
1626 natoms, atoms->atom, ePBC, fr.box, fr.x);
1628 /* Copy the input trxframe struct to the output trxframe struct */
1630 frout.bV = (frout.bV && bVels);
1631 frout.bF = (frout.bF && bForce);
1632 frout.natoms = nout;
1633 if (bNeedPrec && (bSetPrec || !fr.bPrec))
1649 for (i = 0; i < nout; i++)
1651 copy_rvec(fr.x[index[i]], frout.x[i]);
1654 copy_rvec(fr.v[index[i]], frout.v[i]);
1658 copy_rvec(fr.f[index[i]], frout.f[i]);
1663 if (opt2parg_bSet("-shift", NPA, pa))
1665 for (i = 0; i < nout; i++)
1667 for (d = 0; d < DIM; d++)
1669 frout.x[i][d] += outframe*shift[d];
1676 bSplitHere = bSplit && bRmod(fr.time, tzero, split_t);
1680 /* round() is not C89 compatible, so we do this: */
1681 bSplitHere = bSplit && bRmod(floor(fr.time+0.5),
1683 floor(split_t+0.5));
1685 if (bSeparate || bSplitHere)
1687 mk_filenm(outf_base, ftp2ext(ftp), nzero, file_nr, out_file2);
1702 trxout = open_trx(out_file2, filemode);
1709 if (clust_status_id[my_clust] == -1)
1711 sprintf(buf, "%s.%s", clust->grpname[my_clust], ftp2ext(ftp));
1712 clust_status[my_clust] = open_trx(buf, "w");
1713 clust_status_id[my_clust] = 1;
1716 else if (clust_status_id[my_clust] == -2)
1718 gmx_fatal(FARGS, "File %s.xtc should still be open (%d open .xtc files)\n" "in order to write frame %d. my_clust = %d",
1719 clust->grpname[my_clust], ntrxopen, frame,
1722 write_trxframe(clust_status[my_clust], &frout, gc);
1723 nfwritten[my_clust]++;
1724 if (nfwritten[my_clust] ==
1725 (clust->clust->index[my_clust+1]-
1726 clust->clust->index[my_clust]))
1728 close_trx(clust_status[my_clust]);
1729 clust_status[my_clust] = NULL;
1730 clust_status_id[my_clust] = -2;
1734 gmx_fatal(FARGS, "Less than zero open .xtc files!");
1741 write_trxframe(trxout, &frout, gc);
1747 sprintf(title, "Generated by trjconv : %s t= %9.5f",
1748 top_title, fr.time);
1749 if (bSeparate || bSplitHere)
1751 out = ffopen(out_file2, "w");
1756 write_hconf_p(out, title, &useatoms, prec2ndec(frout.prec),
1757 frout.x, frout.bV ? frout.v : NULL, frout.box);
1760 fprintf(out, "REMARK GENERATED BY TRJCONV\n");
1761 sprintf(title, "%s t= %9.5f", top_title, frout.time);
1762 /* if reading from pdb, we want to keep the original
1763 model numbering else we write the output frame
1764 number plus one, because model 0 is not allowed in pdb */
1765 if (ftpin == efPDB && fr.bStep && fr.step > model_nr)
1773 write_pdbfile(out, title, &useatoms, frout.x,
1774 frout.ePBC, frout.box, ' ', model_nr, gc, TRUE);
1777 frout.title = title;
1778 if (bSeparate || bTDump)
1780 frout.bTitle = TRUE;
1783 frout.bAtoms = TRUE;
1785 frout.atoms = &useatoms;
1786 frout.bStep = FALSE;
1787 frout.bTime = FALSE;
1791 frout.bTitle = (outframe == 0);
1792 frout.bAtoms = FALSE;
1796 write_g96_conf(out, &frout, -1, NULL);
1805 gmx_fatal(FARGS, "DHE, ftp=%d\n", ftp);
1807 if (bSeparate || bSplitHere)
1812 /* execute command */
1816 sprintf(c, "%s %d", exec_command, file_nr-1);
1817 /*fprintf(stderr,"Executing '%s'\n",c);*/
1818 #ifdef GMX_NO_SYSTEM
1819 printf("Warning-- No calls to system(3) supported on this platform.");
1820 printf("Warning-- Skipping execution of 'system(\"%s\")'.", c);
1824 gmx_fatal(FARGS, "Error executing command: %s", c);
1832 bHaveNextFrame = read_next_frame(oenv, status, &fr);
1834 while (!(bTDump && bDumpFrame) && bHaveNextFrame);
1837 if (!bHaveFirstFrame || (bTDump && !bDumpFrame))
1839 fprintf(stderr, "\nWARNING no output, "
1840 "last frame read at t=%g\n", fr.time);
1842 fprintf(stderr, "\n");
1849 gmx_rmpbc_done(gpbc);
1856 else if (out != NULL)
1862 for (i = 0; (i < clust->clust->nr); i++)
1864 if (clust_status_id[i] >= 0)
1866 close_trx(clust_status[i]);
1872 do_view(oenv, out_file, NULL);