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
67 enum { euSel,euRect, euTric, euCompact, euNR};
69 static real calc_isquared(int nmol,rvec m_com[],rvec m_shift[],rvec clust_com)
75 for(i=0; (i<nmol); i++) {
76 rvec_add(m_com[i],m_shift[i],m0);
77 rvec_sub(m0,clust_com,dx);
83 static void calc_pbc_cluster(int ecenter,int nrefat,t_topology *top,int ePBC,
84 rvec x[],atom_id index[],
85 rvec clust_com,matrix box)
89 int m,i,j,j0,j1,jj,ai,iter,is;
90 real fac,Isq,min_dist2;
91 rvec dx,ddx,xtest,xrm,box_center;
92 int nmol,nmol_cl,imol_center;
95 rvec *m_com,*m_shift,m0;
98 calc_box_center(ecenter,box,box_center);
100 /* Initiate the pbc structure */
101 memset(&pbc,0,sizeof(pbc));
102 set_pbc(&pbc,ePBC,box);
104 /* Convert atom index to molecular */
106 molind = top->mols.index;
110 snew(bTmp,top->atoms.nr);
112 for(i=0; (i<nrefat); i++) {
113 /* Mark all molecules in the index */
116 /* Binary search assuming the molecules are sorted */
120 if (ai < molind[j0+1])
122 else if (ai >= molind[j1])
126 if (ai < molind[jj+1])
134 /* Double check whether all atoms in all molecules that are marked are part
135 * of the cluster. Simultaneously compute the center of geometry.
137 min_dist2 = 10*sqr(trace(box));
139 for(i=0; (i<nmol); i++) {
140 for(j=molind[i]; (j<molind[i+1]); j++) {
141 if (bMol[i] && !bTmp[j])
142 gmx_fatal(FARGS,"Molecule %d marked for clustering but not atom %d",
144 else if (!bMol[i] && bTmp[j])
145 gmx_fatal(FARGS,"Atom %d marked for clustering but not molecule %d",
148 /* Compute center of geometry of molecule */
149 rvec_inc(m_com[i],x[j]);
153 /* Normalize center of geometry */
154 fac = 1.0/(molind[i+1]-molind[i]);
155 for(m=0; (m<DIM); m++)
157 /* Determine which molecule is closest to the center of the box */
158 pbc_dx(&pbc,box_center,m_com[i],dx);
159 if (iprod(dx,dx) < min_dist2) {
160 min_dist2 = iprod(dx,dx);
169 fprintf(stderr,"No molecules selected in the cluster\n");
171 } else if (imol_center == -1) {
172 fprintf(stderr,"No central molecules could be found\n");
177 /* First calculation is incremental */
178 clear_rvec(clust_com);
180 for(i=m=0; (i<nmol); i++) {
181 /* Check whether this molecule is part of the cluster */
183 if ((i > 0) && (m > 0)) {
184 /* Compute center of cluster by dividing by number of molecules */
185 svmul(1.0/m,clust_com,xrm);
186 /* Distance vector between molecular COM and cluster */
187 pbc_dx(&pbc,m_com[i],xrm,dx);
188 rvec_add(xrm,dx,xtest);
189 /* xtest is now the image of m_com[i] that is closest to clust_com */
190 rvec_inc(clust_com,xtest);
191 rvec_sub(xtest,m_com[i],m_shift[i]);
194 rvec_inc(clust_com,m_com[i]);
199 assert(m == nmol_cl);
200 svmul(1/nmol_cl,clust_com,clust_com);
201 put_atom_in_box(box,clust_com);
203 /* Now check if any molecule is more than half the box from the COM */
208 for(i=0; (i<nmol) && !bChanged; i++) {
210 /* Sum com and shift from com */
211 rvec_add(m_com[i],m_shift[i],m0);
212 pbc_dx(&pbc,m0,clust_com,dx);
213 rvec_add(clust_com,dx,xtest);
214 rvec_sub(xtest,m0,ddx);
215 if (iprod(ddx,ddx) > tol) {
216 /* Here we have used the wrong image for contributing to the COM */
217 rvec_sub(xtest,m_com[i],m_shift[i]);
218 for(j=0; (j<DIM); j++)
219 clust_com[j] += (xtest[j]-m0[j])/nmol_cl;
224 Isq = calc_isquared(nmol,m_com,m_shift,clust_com);
225 put_atom_in_box(box,clust_com);
227 if (bChanged && (iter > 0))
228 printf("COM: %8.3f %8.3f %8.3f iter = %d Isq = %8.3f\n",
229 clust_com[XX],clust_com[YY],clust_com[ZZ],iter,Isq);
233 /* Now transfer the shift to all atoms in the molecule */
234 for(i=0; (i<nmol); i++) {
236 for(j=molind[i]; (j<molind[i+1]); j++)
237 rvec_inc(x[j],m_shift[i]);
245 static void put_molecule_com_in_box(int unitcell_enum,int ecenter,
247 int natoms,t_atom atom[],
248 int ePBC,matrix box,rvec x[])
252 rvec com,new_com,shift,dx,box_center;
257 calc_box_center(ecenter,box,box_center);
258 set_pbc(&pbc,ePBC,box);
260 gmx_fatal(FARGS,"There are no molecule descriptions. I need a tpr file for this pbc option.");
261 for(i=0; (i<mols->nr); i++) {
265 for(j=mols->index[i]; (j<mols->index[i+1] && j<natoms); j++) {
271 /* calculate final COM */
272 svmul(1.0/mtot, com, com);
274 /* check if COM is outside box */
275 copy_rvec(com,new_com);
276 switch (unitcell_enum) {
278 put_atoms_in_box(box,1,&new_com);
281 put_atoms_in_triclinic_unitcell(ecenter,box,1,&new_com);
284 put_atoms_in_compact_unitcell(ePBC,ecenter,box,1,&new_com);
287 rvec_sub(new_com,com,shift);
288 if (norm2(shift) > 0) {
290 fprintf (debug,"\nShifting position of molecule %d "
291 "by %8.3f %8.3f %8.3f\n", i+1, PR_VEC(shift));
292 for(j=mols->index[i]; (j<mols->index[i+1] && j<natoms); j++) {
293 rvec_inc(x[j],shift);
299 static void put_residue_com_in_box(int unitcell_enum,int ecenter,
300 int natoms,t_atom atom[],
301 int ePBC,matrix box,rvec x[])
303 atom_id i, j, res_start, res_end, res_nat;
307 rvec box_center,com,new_com,shift;
309 calc_box_center(ecenter,box,box_center);
315 for(i=0; i<natoms+1; i++) {
316 if (i == natoms || (presnr != atom[i].resind && presnr != NOTSET)) {
317 /* calculate final COM */
319 res_nat = res_end - res_start;
320 svmul(1.0/mtot, com, com);
322 /* check if COM is outside box */
323 copy_rvec(com,new_com);
324 switch (unitcell_enum) {
326 put_atoms_in_box(box,1,&new_com);
329 put_atoms_in_triclinic_unitcell(ecenter,box,1,&new_com);
332 put_atoms_in_compact_unitcell(ePBC,ecenter,box,1,&new_com);
335 rvec_sub(new_com,com,shift);
338 fprintf (debug,"\nShifting position of residue %d (atoms %u-%u) "
339 "by %g,%g,%g\n", atom[res_start].resind+1,
340 res_start+1, res_end+1, PR_VEC(shift));
341 for(j=res_start; j<res_end; j++)
342 rvec_inc(x[j],shift);
347 /* remember start of new residue */
357 presnr = atom[i].resind;
362 static void center_x(int ecenter,rvec x[],matrix box,int n,int nc,atom_id ci[])
365 rvec cmin,cmax,box_center,dx;
368 copy_rvec(x[ci[0]],cmin);
369 copy_rvec(x[ci[0]],cmax);
370 for(i=0; i<nc; i++) {
372 for(m=0; m<DIM; m++) {
373 if (x[ai][m] < cmin[m])
375 else if (x[ai][m] > cmax[m])
379 calc_box_center(ecenter,box,box_center);
381 dx[m] = box_center[m]-(cmin[m]+cmax[m])*0.5;
388 static void mk_filenm(char *base,const char *ext,int ndigit,int file_nr,
394 strcpy(out_file,base);
401 strncat(out_file,"00000000000",ndigit-nd);
402 sprintf(nbuf,"%d.",file_nr);
403 strcat(out_file,nbuf);
404 strcat(out_file,ext);
407 void check_trn(const char *fn)
409 if ((fn2ftp(fn) != efTRJ) && (fn2ftp(fn) != efTRR))
410 gmx_fatal(FARGS,"%s is not a trj file, exiting\n",fn);
413 #if (!defined WIN32 && !defined _WIN32 && !defined WIN64 && !defined _WIN64)
414 void do_trunc(const char *fn, real t0)
426 gmx_fatal(FARGS,"You forgot to set the truncation time");
428 /* Check whether this is a .trj file */
431 in = open_trn(fn,"r");
432 fp = gmx_fio_getfp(in);
434 fprintf(stderr,"Sorry, can not trunc %s, truncation of this filetype is not supported\n",fn);
438 fpos = gmx_fio_ftell(in);
440 while (!bStop && fread_trnheader(in,&sh,&bOK)) {
441 fread_htrn(in,&sh,NULL,NULL,NULL,NULL);
445 gmx_fseek(fp, fpos, SEEK_SET);
450 fprintf(stderr,"Do you REALLY want to truncate this trajectory (%s) at:\n"
451 "frame %d, time %g, bytes %ld ??? (type YES if so)\n",
452 fn,j,t,(long int)fpos);
453 if(1 != scanf("%s",yesno))
455 gmx_fatal(FARGS,"Error reading user input");
457 if (strcmp(yesno,"YES") == 0) {
458 fprintf(stderr,"Once again, I'm gonna DO this...\n");
460 if(0 != truncate(fn,fpos))
462 gmx_fatal(FARGS,"Error truncating file %s",fn);
466 fprintf(stderr,"Ok, I'll forget about it\n");
470 fprintf(stderr,"Already at end of file (t=%g)...\n",t);
477 int gmx_trjconv(int argc,char *argv[])
479 const char *desc[] = {
480 "trjconv can convert trajectory files in many ways:[BR]",
481 "[BB]1.[bb] from one format to another[BR]",
482 "[BB]2.[bb] select a subset of atoms[BR]",
483 "[BB]3.[bb] change the periodicity representation[BR]",
484 "[BB]4.[bb] keep multimeric molecules together[BR]",
485 "[BB]5.[bb] center atoms in the box[BR]",
486 "[BB]6.[bb] fit atoms to reference structure[BR]",
487 "[BB]7.[bb] reduce the number of frames[BR]",
488 "[BB]8.[bb] change the timestamps of the frames ",
489 "([TT]-t0[tt] and [TT]-timestep[tt])[BR]",
490 "[BB]9.[bb] cut the trajectory in small subtrajectories according",
491 "to information in an index file. This allows subsequent analysis of",
492 "the subtrajectories that could, for example be the result of a",
493 "cluster analysis. Use option [TT]-sub[tt].",
494 "This assumes that the entries in the index file are frame numbers and",
495 "dumps each group in the index file to a separate trajectory file.[BR]",
496 "[BB]10.[bb] select frames within a certain range of a quantity given",
497 "in an [TT].xvg[tt] file.[PAR]",
498 "The program [TT]trjcat[tt] can concatenate multiple trajectory files.",
500 "Currently seven formats are supported for input and output:",
501 "[TT].xtc[tt], [TT].trr[tt], [TT].trj[tt], [TT].gro[tt], [TT].g96[tt],",
502 "[TT].pdb[tt] and [TT].g87[tt].",
503 "The file formats are detected from the file extension.",
504 "The precision of [TT].xtc[tt] and [TT].gro[tt] output is taken from the",
505 "input file for [TT].xtc[tt], [TT].gro[tt] and [TT].pdb[tt],",
506 "and from the [TT]-ndec[tt] option for other input formats. The precision",
507 "is always taken from [TT]-ndec[tt], when this option is set.",
508 "All other formats have fixed precision. [TT].trr[tt] and [TT].trj[tt]",
509 "output can be single or double precision, depending on the precision",
510 "of the trjconv binary.",
511 "Note that velocities are only supported in",
512 "[TT].trr[tt], [TT].trj[tt], [TT].gro[tt] and [TT].g96[tt] files.[PAR]",
513 "Option [TT]-app[tt] can be used to",
514 "append output to an existing trajectory file.",
515 "No checks are performed to ensure integrity",
516 "of the resulting combined trajectory file.[PAR]",
517 "Option [TT]-sep[tt] can be used to write every frame to a separate",
518 ".gro, .g96 or .pdb file, default all frames all written to one file.",
519 "[TT].pdb[tt] files with all frames concatenated can be viewed with",
520 "[TT]rasmol -nmrpdb[tt].[PAR]",
521 "It is possible to select part of your trajectory and write it out",
522 "to a new trajectory file in order to save disk space, e.g. for leaving",
523 "out the water from a trajectory of a protein in water.",
524 "[BB]ALWAYS[bb] put the original trajectory on tape!",
525 "We recommend to use the portable [TT].xtc[tt] format for your analysis",
526 "to save disk space and to have portable files.[PAR]",
527 "There are two options for fitting the trajectory to a reference",
528 "either for essential dynamics analysis or for whatever.",
529 "The first option is just plain fitting to a reference structure",
530 "in the structure file, the second option is a progressive fit",
531 "in which the first timeframe is fitted to the reference structure ",
532 "in the structure file to obtain and each subsequent timeframe is ",
533 "fitted to the previously fitted structure. This way a continuous",
534 "trajectory is generated, which might not be the case when using the",
535 "regular fit method, e.g. when your protein undergoes large",
536 "conformational transitions.[PAR]",
537 "Option [TT]-pbc[tt] sets the type of periodic boundary condition",
539 "* [TT]mol[tt] puts the center of mass of molecules in the box.[BR]",
540 "* [TT]res[tt] puts the center of mass of residues in the box.[BR]",
541 "* [TT]atom[tt] puts all the atoms in the box.[BR]",
542 "* [TT]nojump[tt] checks if atoms jump across the box and then puts",
543 "them back. This has the effect that all molecules",
544 "will remain whole (provided they were whole in the initial",
545 "conformation), note that this ensures a continuous trajectory but",
546 "molecules may diffuse out of the box. The starting configuration",
547 "for this procedure is taken from the structure file, if one is",
548 "supplied, otherwise it is the first frame.[BR]",
549 "* [TT]cluster[tt] clusters all the atoms in the selected index",
550 "such that they are all closest to the center of mass of the cluster",
551 "which is iteratively updated. Note that this will only give meaningful",
552 "results if you in fact have a cluster. Luckily that can be checked",
553 "afterwards using a trajectory viewer. Note also that if your molecules",
554 "are broken this will not work either.[BR]",
555 "* [TT]whole[tt] only makes broken molecules whole.[PAR]",
556 "Option [TT]-ur[tt] sets the unit cell representation for options",
557 "[TT]mol[tt], [TT]res[tt] and [TT]atom[tt] of [TT]-pbc[tt].",
558 "All three options give different results for triclinic boxes and",
559 "identical results for rectangular boxes.",
560 "[TT]rect[tt] is the ordinary brick shape.",
561 "[TT]tric[tt] is the triclinic unit cell.",
562 "[TT]compact[tt] puts all atoms at the closest distance from the center",
563 "of the box. This can be useful for visualizing e.g. truncated",
564 "octahedrons. The center for options [TT]tric[tt] and [TT]compact[tt]",
565 "is [TT]tric[tt] (see below), unless the option [TT]-boxcenter[tt]",
566 "is set differently.[PAR]",
567 "Option [TT]-center[tt] centers the system in the box. The user can",
568 "select the group which is used to determine the geometrical center.",
569 "Option [TT]-boxcenter[tt] sets the location of the center of the box",
570 "for options [TT]-pbc[tt] and [TT]-center[tt]. The center options are:",
571 "[TT]tric[tt]: half of the sum of the box vectors,",
572 "[TT]rect[tt]: half of the box diagonal,",
573 "[TT]zero[tt]: zero.",
574 "Use option [TT]-pbc mol[tt] in addition to [TT]-center[tt] when you",
575 "want all molecules in the box after the centering.[PAR]",
576 "With [TT]-dt[tt] it is possible to reduce the number of ",
577 "frames in the output. This option relies on the accuracy of the times",
578 "in your input trajectory, so if these are inaccurate use the",
579 "[TT]-timestep[tt] option to modify the time (this can be done",
580 "simultaneously). For making smooth movies the program [TT]g_filter[tt]",
581 "can reduce the number of frames while using low-pass frequency",
582 "filtering, this reduces aliasing of high frequency motions.[PAR]",
583 "Using [TT]-trunc[tt] trjconv can truncate [TT].trj[tt] in place, i.e.",
584 "without copying the file. This is useful when a run has crashed",
585 "during disk I/O (one more disk full), or when two contiguous",
586 "trajectories must be concatenated without have double frames.[PAR]",
587 "[TT]trjcat[tt] is more suitable for concatenating trajectory files.[PAR]",
588 "Option [TT]-dump[tt] can be used to extract a frame at or near",
589 "one specific time from your trajectory.[PAR]",
590 "Option [TT]-drop[tt] reads an [TT].xvg[tt] file with times and values.",
591 "When options [TT]-dropunder[tt] and/or [TT]-dropover[tt] are set,",
592 "frames with a value below and above the value of the respective options",
593 "will not be written."
609 const char *pbc_opt[epNR + 1] =
610 { NULL, "none", "mol", "res", "atom", "nojump", "cluster", "whole",
614 const char *unitcell_opt[euNR+1] =
615 { NULL, "rect", "tric", "compact", NULL };
618 { ecSel, ecTric, ecRect, ecZero, ecNR};
619 const char *center_opt[ecNR+1] =
620 { NULL, "tric", "rect", "zero", NULL };
626 efSel, efNone, efFit, efFitXY, efReset, efResetXY, efPFit, efNR
628 const char *fit[efNR + 1] =
629 { NULL, "none", "rot+trans", "rotxy+transxy", "translation", "transxy",
630 "progressive", NULL };
632 static gmx_bool bAppend=FALSE,bSeparate=FALSE,bVels=TRUE,bForce=FALSE,bCONECT=FALSE;
633 static gmx_bool bCenter=FALSE;
634 static int skip_nr=1,ndec=3,nzero=0;
635 static real tzero=0,delta_t=0,timestep=0,ttrunc=-1,tdump=-1,split_t=0;
636 static rvec newbox = {0,0,0}, shift = {0,0,0}, trans = {0,0,0};
637 static char *exec_command=NULL;
638 static real dropunder=0,dropover=0;
639 static gmx_bool bRound=FALSE;
644 { "-skip", FALSE, etINT,
645 { &skip_nr }, "Only write every nr-th frame" },
646 { "-dt", FALSE, etTIME,
648 "Only write frame when t MOD dt = first time (%t)" },
649 { "-round", FALSE, etBOOL,
650 { &bRound }, "Round measurements to nearest picosecond"
652 { "-dump", FALSE, etTIME,
653 { &tdump }, "Dump frame nearest specified time (%t)" },
654 { "-t0", FALSE, etTIME,
656 "Starting time (%t) (default: don't change)" },
657 { "-timestep", FALSE, etTIME,
659 "Change time step between input frames (%t)" },
660 { "-pbc", FALSE, etENUM,
662 "PBC treatment (see help text for full description)" },
663 { "-ur", FALSE, etENUM,
664 { unitcell_opt }, "Unit-cell representation" },
665 { "-center", FALSE, etBOOL,
666 { &bCenter }, "Center atoms in box" },
667 { "-boxcenter", FALSE, etENUM,
668 { center_opt }, "Center for -pbc and -center" },
669 { "-box", FALSE, etRVEC,
671 "Size for new cubic box (default: read from input)" },
672 { "-trans", FALSE, etRVEC,
674 "All coordinates will be translated by trans. This "
675 "can advantageously be combined with -pbc mol -ur "
677 { "-shift", FALSE, etRVEC,
679 "All coordinates will be shifted by framenr*shift" },
680 { "-fit", FALSE, etENUM,
682 "Fit molecule to ref structure in the structure file" },
683 { "-ndec", FALSE, etINT,
685 "Precision for .xtc and .gro writing in number of "
687 { "-vel", FALSE, etBOOL,
688 { &bVels }, "Read and write velocities if possible" },
689 { "-force", FALSE, etBOOL,
690 { &bForce }, "Read and write forces if possible" },
691 #if (!defined WIN32 && !defined _WIN32 && !defined WIN64 && !defined _WIN64)
692 { "-trunc", FALSE, etTIME,
694 "Truncate input trj file after this time (%t)" },
696 { "-exec", FALSE, etSTR,
698 "Execute command for every output frame with the "
699 "frame number as argument" },
700 { "-app", FALSE, etBOOL,
701 { &bAppend }, "Append output" },
702 { "-split", FALSE, etTIME,
704 "Start writing new file when t MOD split = first "
706 { "-sep", FALSE, etBOOL,
708 "Write each frame to a separate .gro, .g96 or .pdb "
710 { "-nzero", FALSE, etINT,
712 "Prepend file number in case you use the -sep flag "
713 "with this number of zeroes" },
714 { "-dropunder", FALSE, etREAL,
715 { &dropunder }, "Drop all frames below this value" },
716 { "-dropover", FALSE, etREAL,
717 { &dropover }, "Drop all frames above this value" },
718 { "-conect", FALSE, etBOOL,
720 "Add conect records when writing pdb files. Useful "
721 "for visualization of non-standard molecules, e.g. "
722 "coarse grained ones" } };
723 #define NPA asize(pa)
726 t_trxstatus *trxout=NULL;
728 int ftp,ftpin=0,file_nr;
731 rvec *xmem=NULL,*vmem=NULL,*fmem=NULL;
732 rvec *xp=NULL,x_shift,hbox,box_center,dx;
733 real xtcpr, lambda,*w_rls=NULL;
734 int m,i,d,frame,outframe,natoms,nout,ncent,nre,newstep=0,model_nr;
739 t_atoms *atoms=NULL,useatoms;
741 atom_id *index,*cindex;
745 int ifit,irms,my_clust=-1;
746 atom_id *ind_fit,*ind_rms;
747 char *gn_fit,*gn_rms;
748 t_cluster_ndx *clust=NULL;
749 t_trxstatus **clust_status=NULL;
750 int *clust_status_id=NULL;
753 int ndrop=0,ncol,drop0=0,drop1=0,dropuse=0;
755 real tshift=0,t0=-1,dt=0.001,prec;
756 gmx_bool bFit,bFitXY,bPFit,bReset;
758 gmx_rmpbc_t gpbc=NULL;
759 gmx_bool bRmPBC,bPBCWhole,bPBCcomRes,bPBCcomMol,bPBCcomAtom,bPBC,bNoJump,bCluster;
760 gmx_bool bCopy,bDoIt,bIndex,bTDump,bSetTime,bTPS=FALSE,bDTset=FALSE;
761 gmx_bool bExec,bTimeStep=FALSE,bDumpFrame=FALSE,bSetPrec,bNeedPrec;
762 gmx_bool bHaveFirstFrame,bHaveNextFrame,bSetBox,bSetUR,bSplit=FALSE;
763 gmx_bool bSubTraj=FALSE,bDropUnder=FALSE,bDropOver=FALSE,bTrans=FALSE;
764 gmx_bool bWriteFrame,bSplitHere;
765 const char *top_file,*in_file,*out_file=NULL;
766 char out_file2[256],*charpt;
767 char *outf_base=NULL;
768 const char *outf_ext=NULL;
769 char top_title[256],title[256],command[256],filemode[5];
771 gmx_bool bWarnCompact=FALSE;
776 { efTRX, "-f", NULL, ffREAD },
777 { efTRO, "-o", NULL, ffWRITE },
778 { efTPS, NULL, NULL, ffOPTRD },
779 { efNDX, NULL, NULL, ffOPTRD },
780 { efNDX, "-fr", "frames", ffOPTRD },
781 { efNDX, "-sub", "cluster", ffOPTRD },
782 { efXVG, "-drop","drop", ffOPTRD }
784 #define NFILE asize(fnm)
786 CopyRight(stderr,argv[0]);
787 parse_common_args(&argc,argv,
788 PCA_CAN_BEGIN | PCA_CAN_END | PCA_CAN_VIEW |
789 PCA_TIME_UNIT | PCA_BE_NICE,
790 NFILE,fnm,NPA,pa,asize(desc),desc,
793 top_file = ftp2fn(efTPS,NFILE,fnm);
796 /* Check command line */
797 in_file=opt2fn("-f",NFILE,fnm);
800 #if (!defined WIN32 && !defined _WIN32 && !defined WIN64 && !defined _WIN64)
801 do_trunc(in_file,ttrunc);
805 /* mark active cmdline options */
806 bSetBox = opt2parg_bSet("-box", NPA, pa);
807 bSetTime = opt2parg_bSet("-t0", NPA, pa);
808 bSetPrec = opt2parg_bSet("-ndec", NPA, pa);
809 bSetUR = opt2parg_bSet("-ur", NPA, pa);
810 bExec = opt2parg_bSet("-exec", NPA, pa);
811 bTimeStep = opt2parg_bSet("-timestep", NPA, pa);
812 bTDump = opt2parg_bSet("-dump", NPA, pa);
813 bDropUnder= opt2parg_bSet("-dropunder", NPA, pa);
814 bDropOver = opt2parg_bSet("-dropover", NPA, pa);
815 bTrans = opt2parg_bSet("-trans",NPA,pa);
816 bSplit = (split_t != 0);
818 /* parse enum options */
819 fit_enum = nenum(fit);
820 bFit = (fit_enum==efFit || fit_enum==efFitXY);
821 bFitXY = fit_enum==efFitXY;
822 bReset = (fit_enum==efReset || fit_enum==efResetXY);
823 bPFit = fit_enum==efPFit;
824 pbc_enum = nenum(pbc_opt);
825 bPBCWhole = pbc_enum==epWhole;
826 bPBCcomRes = pbc_enum==epComRes;
827 bPBCcomMol = pbc_enum==epComMol;
828 bPBCcomAtom= pbc_enum==epComAtom;
829 bNoJump = pbc_enum==epNojump;
830 bCluster = pbc_enum==epCluster;
831 bPBC = pbc_enum!=epNone;
832 unitcell_enum = nenum(unitcell_opt);
833 ecenter = nenum(center_opt) - ecTric;
835 /* set and check option dependencies */
836 if (bPFit) bFit = TRUE; /* for pfit, fit *must* be set */
837 if (bFit) bReset = TRUE; /* for fit, reset *must* be set */
840 nfitdim = (fit_enum==efFitXY || fit_enum==efResetXY) ? 2 : 3;
841 bRmPBC = bFit || bPBCWhole || bPBCcomRes || bPBCcomMol;
844 if (!(bPBCcomRes || bPBCcomMol || bPBCcomAtom)) {
846 "WARNING: Option for unitcell representation (-ur %s)\n"
847 " only has effect in combination with -pbc %s, %s or %s.\n"
848 " Ingoring unitcell representation.\n\n",
849 unitcell_opt[0],pbc_opt[2],pbc_opt[3],pbc_opt[4]);
854 gmx_fatal(FARGS,"PBC condition treatment does not work together with rotational fit.\n"
855 "Please do the PBC condition treatment first and then run trjconv in a second step\n"
856 "for the rotational fit.\n"
857 "First doing the rotational fit and then doing the PBC treatment gives incorrect\n"
861 /* ndec is in nr of decimal places, prec is a multiplication factor: */
863 for (i=0; i<ndec; i++)
866 bIndex=ftp2bSet(efNDX,NFILE,fnm);
869 /* Determine output type */
870 out_file=opt2fn("-o",NFILE,fnm);
871 ftp=fn2ftp(out_file);
872 fprintf(stderr,"Will write %s: %s\n",ftp2ext(ftp),ftp2desc(ftp));
873 bNeedPrec = (ftp==efXTC || ftp==efGRO);
875 /* check if velocities are possible in input and output files */
876 ftpin=fn2ftp(in_file);
877 bVels= (ftp==efTRR || ftp==efTRJ || ftp==efGRO || ftp==efG96)
878 && (ftpin==efTRR || ftpin==efTRJ || ftpin==efGRO || ftpin==efG96 ||
881 if (bSeparate || bSplit) {
882 outf_ext = strrchr(out_file,'.');
883 if (outf_ext == NULL)
884 gmx_fatal(FARGS,"Output file name '%s' does not contain a '.'",out_file);
885 outf_base = strdup(out_file);
886 outf_base[outf_ext - out_file] = '\0';
889 bSubTraj = opt2bSet("-sub",NFILE,fnm);
891 if ((ftp != efXTC) && (ftp != efTRN))
892 gmx_fatal(FARGS,"Can only use the sub option with output file types "
894 clust = cluster_index(NULL,opt2fn("-sub",NFILE,fnm));
896 /* Check for number of files disabled, as FOPEN_MAX is not the correct
897 * number to check for. In my linux box it is only 16.
899 if (0 && (clust->clust->nr > FOPEN_MAX-4))
900 gmx_fatal(FARGS,"Can not open enough (%d) files to write all the"
901 " trajectories.\ntry splitting the index file in %d parts.\n"
903 clust->clust->nr,1+clust->clust->nr/FOPEN_MAX,FOPEN_MAX);
905 snew(clust_status,clust->clust->nr);
906 snew(clust_status_id,clust->clust->nr);
907 snew(nfwritten,clust->clust->nr);
908 for(i=0; (i<clust->clust->nr); i++)
910 clust_status[i] = NULL;
911 clust_status_id[i] = -1;
913 bSeparate = bSplit = FALSE;
919 /* Determine whether to read a topology */
920 bTPS = (ftp2bSet(efTPS,NFILE,fnm) ||
921 bRmPBC || bReset || bPBCcomMol ||
922 (ftp == efGRO) || (ftp == efPDB) || bCONECT);
924 /* Determine if when can read index groups */
925 bIndex = (bIndex || bTPS);
928 read_tps_conf(top_file,top_title,&top,&ePBC,&xp,NULL,top_box,
929 bReset || bPBCcomRes);
931 /* top_title is only used for gro and pdb,
932 * the header in such a file is top_title t= ...
933 * to prevent a double t=, remove it from top_title
935 if ((charpt=strstr(top_title," t= ")))
939 gc = gmx_conect_generate(&top);
941 gpbc = gmx_rmpbc_init(&top.idef,ePBC,top.atoms.nr,top_box);
944 /* get frame number index */
946 if (opt2bSet("-fr",NFILE,fnm)) {
947 printf("Select groups of frame number indices:\n");
948 rd_index(opt2fn("-fr",NFILE,fnm),1,&nrfri,(atom_id **)&frindex,&frname);
950 for(i=0; i<nrfri; i++)
951 fprintf(debug,"frindex[%4d]=%4d\n",i,frindex[i]);
954 /* get index groups etc. */
956 printf("Select group for %s fit\n",
957 bFit?"least squares":"translational");
958 get_index(atoms,ftp2fn_null(efNDX,NFILE,fnm),
959 1,&ifit,&ind_fit,&gn_fit);
963 gmx_fatal(FARGS,"Need at least 2 atoms to fit!\n");
965 fprintf(stderr,"WARNING: fitting with only 2 atoms is not unique\n");
969 printf("Select group for clustering\n");
970 get_index(atoms,ftp2fn_null(efNDX,NFILE,fnm),
971 1,&ifit,&ind_fit,&gn_fit);
976 printf("Select group for centering\n");
977 get_index(atoms,ftp2fn_null(efNDX,NFILE,fnm),
978 1,&ncent,&cindex,&grpnm);
980 printf("Select group for output\n");
981 get_index(atoms,ftp2fn_null(efNDX,NFILE,fnm),
982 1,&nout,&index,&grpnm);
984 /* no index file, so read natoms from TRX */
985 if (!read_first_frame(oenv,&status,in_file,&fr,TRX_DONT_SKIP))
986 gmx_fatal(FARGS,"Could not read a frame from %s",in_file);
991 for(i=0;i<natoms;i++)
1001 snew(w_rls,atoms->nr);
1002 for(i=0; (i<ifit); i++)
1003 w_rls[ind_fit[i]]=atoms->atom[ind_fit[i]].m;
1005 /* Restore reference structure and set to origin,
1006 store original location (to put structure back) */
1008 gmx_rmpbc(gpbc,top.atoms.nr,top_box,xp);
1009 copy_rvec(xp[index[0]],x_shift);
1010 reset_x_ndim(nfitdim,ifit,ind_fit,atoms->nr,NULL,xp,w_rls);
1011 rvec_dec(x_shift,xp[index[0]]);
1013 clear_rvec(x_shift);
1015 if (bDropUnder || bDropOver) {
1016 /* Read the xvg file with the drop values */
1017 fprintf(stderr,"\nReading drop file ...");
1018 ndrop = read_xvg(opt2fn("-drop",NFILE,fnm),&dropval,&ncol);
1019 fprintf(stderr," %d time points\n",ndrop);
1020 if (ndrop == 0 || ncol < 2)
1021 gmx_fatal(FARGS,"Found no data points in %s",
1022 opt2fn("-drop",NFILE,fnm));
1027 /* Make atoms struct for output in GRO or PDB files */
1028 if ((ftp == efGRO) || ((ftp == efG96) && bTPS) || (ftp == efPDB)) {
1029 /* get memory for stuff to go in pdb file */
1030 init_t_atoms(&useatoms,atoms->nr,FALSE);
1031 sfree(useatoms.resinfo);
1032 useatoms.resinfo = atoms->resinfo;
1033 for(i=0;(i<nout);i++) {
1034 useatoms.atomname[i]=atoms->atomname[index[i]];
1035 useatoms.atom[i]=atoms->atom[index[i]];
1036 useatoms.nres=max(useatoms.nres,useatoms.atom[i].resind+1);
1040 /* select what to read */
1041 if (ftp==efTRR || ftp==efTRJ)
1046 flags = flags | TRX_READ_V;
1048 flags = flags | TRX_READ_F;
1050 /* open trx file for reading */
1051 bHaveFirstFrame = read_first_frame(oenv,&status,in_file,&fr,flags);
1053 fprintf(stderr,"\nPrecision of %s is %g (nm)\n",in_file,1/fr.prec);
1055 if (bSetPrec || !fr.bPrec)
1056 fprintf(stderr,"\nSetting output precision to %g (nm)\n",1/prec);
1058 fprintf(stderr,"Using output precision of %g (nm)\n",1/prec);
1061 if (bHaveFirstFrame) {
1062 set_trxframe_ePBC(&fr,ePBC);
1067 tshift=tzero-fr.time;
1071 /* open output for writing */
1072 if ((bAppend) && (gmx_fexist(out_file))) {
1073 strcpy(filemode,"a");
1074 fprintf(stderr,"APPENDING to existing file %s\n",out_file);
1076 strcpy(filemode,"w");
1083 if (!bSplit && !bSubTraj)
1084 trxout = open_trx(out_file,filemode);
1089 if (( !bSeparate && !bSplit ) && !bSubTraj)
1090 out=ffopen(out_file,filemode);
1096 /* check if index is meaningful */
1097 for(i=0; i<nout; i++) {
1098 if (index[i] >= natoms)
1099 gmx_fatal(FARGS,"Index[%d] %d is larger than the number of atoms in the trajectory file (%d)",i,index[i]+1,natoms);
1100 bCopy = bCopy || (i != index[i]);
1113 fprintf(gmx_fio_getfp(trx_get_fileio(trxout)),
1114 "Generated by %s. #atoms=%d, a BOX is"
1115 " stored in this file.\n",Program(),nout);
1117 /* Start the big loop over frames */
1124 /* Main loop over frames */
1132 /*if (frame >= clust->clust->nra)
1133 gmx_fatal(FARGS,"There are more frames in the trajectory than in the cluster index file\n");*/
1134 if (frame >= clust->maxframe)
1137 my_clust = clust->inv_clust[frame];
1138 if ((my_clust < 0) || (my_clust >= clust->clust->nr) ||
1139 (my_clust == NO_ATID))
1144 /* generate new box */
1146 for (m=0; m<DIM; m++)
1147 fr.box[m][m] = newbox[m];
1151 for(i=0; i<natoms; i++)
1152 rvec_inc(fr.x[i],trans);
1156 /* determine timestep */
1165 /* This is not very elegant, as one can not dump a frame after
1166 * a timestep with is more than twice as small as the first one. */
1167 bDumpFrame = (fr.time > tdump-0.5*dt) && (fr.time <= tdump+0.5*dt);
1171 /* determine if an atom jumped across the box and reset it if so */
1172 if (bNoJump && (bTPS || frame!=0)) {
1173 for(d=0; d<DIM; d++)
1174 hbox[d] = 0.5*fr.box[d][d];
1175 for(i=0; i<natoms; i++) {
1177 rvec_dec(fr.x[i],x_shift);
1178 for(m=DIM-1; m>=0; m--)
1180 while (fr.x[i][m]-xp[i][m] <= -hbox[m])
1182 fr.x[i][d] += fr.box[m][d];
1183 while (fr.x[i][m]-xp[i][m] > hbox[m])
1185 fr.x[i][d] -= fr.box[m][d];
1189 else if (bCluster && bTPS) {
1192 calc_pbc_cluster(ecenter,ifit,&top,ePBC,fr.x,ind_fit,com,fr.box);
1196 /* Now modify the coords according to the flags,
1197 for normal fit, this is only done for output frames */
1199 gmx_rmpbc_trxfr(gpbc,&fr);
1201 reset_x_ndim(nfitdim,ifit,ind_fit,natoms,NULL,fr.x,w_rls);
1202 do_fit(natoms,w_rls,xp,fr.x);
1205 /* store this set of coordinates for future use */
1206 if (bPFit || bNoJump) {
1209 for(i=0; (i<natoms); i++) {
1210 copy_rvec(fr.x[i],xp[i]);
1211 rvec_inc(fr.x[i],x_shift);
1216 /* see if we have a frame from the frame index group */
1217 for(i=0; i<nrfri && !bDumpFrame; i++)
1218 bDumpFrame = frame == frindex[i];
1219 if (debug && bDumpFrame)
1220 fprintf(debug,"dumping %d\n",frame);
1223 ( ( !bTDump && !frindex && frame % skip_nr == 0 ) || bDumpFrame );
1225 if (bWriteFrame && (bDropUnder || bDropOver)) {
1226 while (dropval[0][drop1]<fr.time && drop1+1<ndrop) {
1230 if (fabs(dropval[0][drop0] - fr.time)
1231 < fabs(dropval[0][drop1] - fr.time)) {
1236 if ((bDropUnder && dropval[1][dropuse] < dropunder) ||
1237 (bDropOver && dropval[1][dropuse] > dropover))
1238 bWriteFrame = FALSE;
1245 fr.time = tzero+frame*timestep;
1251 fprintf(stderr,"\nDumping frame at t= %g %s\n",
1252 output_env_conv_time(oenv,fr.time),output_env_get_time_unit(oenv));
1254 /* check for writing at each delta_t */
1255 bDoIt=(delta_t == 0);
1259 bDoIt=bRmod(fr.time,tzero, delta_t);
1261 /* round() is not C89 compatible, so we do this: */
1262 bDoIt=bRmod(floor(fr.time+0.5),floor(tzero+0.5),
1263 floor(delta_t+0.5));
1266 if (bDoIt || bTDump) {
1267 /* print sometimes */
1268 if ( ((outframe % SKIP) == 0) || (outframe < SKIP) )
1269 fprintf(stderr," -> frame %6d time %8.3f \r",
1270 outframe,output_env_conv_time(oenv,fr.time));
1273 /* Now modify the coords according to the flags,
1274 for PFit we did this already! */
1277 gmx_rmpbc_trxfr(gpbc,&fr);
1280 reset_x_ndim(nfitdim,ifit,ind_fit,natoms,NULL,fr.x,w_rls);
1282 do_fit_ndim(nfitdim,natoms,w_rls,xp,fr.x);
1284 for(i=0; i<natoms; i++)
1285 rvec_inc(fr.x[i],x_shift);
1289 center_x(ecenter,fr.x,fr.box,natoms,ncent,cindex);
1293 switch (unitcell_enum) {
1295 put_atoms_in_box(fr.box,natoms,fr.x);
1298 put_atoms_in_triclinic_unitcell(ecenter,fr.box,natoms,fr.x);
1301 warn = put_atoms_in_compact_unitcell(ePBC,ecenter,fr.box,
1303 if (warn && !bWarnCompact) {
1304 fprintf(stderr,"\n%s\n",warn);
1305 bWarnCompact = TRUE;
1311 put_residue_com_in_box(unitcell_enum,ecenter,
1312 natoms,atoms->atom,ePBC,fr.box,fr.x);
1315 put_molecule_com_in_box(unitcell_enum,ecenter,
1317 natoms,atoms->atom,ePBC,fr.box,fr.x);
1319 /* Copy the input trxframe struct to the output trxframe struct */
1321 frout.natoms = nout;
1322 if (bNeedPrec && (bSetPrec || !fr.bPrec)) {
1334 for(i=0; i<nout; i++) {
1335 copy_rvec(fr.x[index[i]],frout.x[i]);
1336 if (bVels && fr.bV) {
1337 copy_rvec(fr.v[index[i]],frout.v[i]);
1339 if (bForce && fr.bF) {
1340 copy_rvec(fr.f[index[i]],frout.f[i]);
1345 if (opt2parg_bSet("-shift",NPA,pa))
1346 for(i=0; i<nout; i++)
1347 for (d=0; d<DIM; d++)
1348 frout.x[i][d] += outframe*shift[d];
1351 bSplitHere = bSplit && bRmod(fr.time,tzero, split_t);
1354 /* round() is not C89 compatible, so we do this: */
1355 bSplitHere = bSplit && bRmod(floor(fr.time+0.5),
1357 floor(split_t+0.5));
1359 if (bSeparate || bSplitHere)
1360 mk_filenm(outf_base,ftp2ext(ftp),nzero,file_nr,out_file2);
1370 trxout = open_trx(out_file2,filemode);
1373 if (my_clust != -1) {
1375 if (clust_status_id[my_clust] == -1) {
1376 sprintf(buf,"%s.%s",clust->grpname[my_clust],ftp2ext(ftp));
1377 clust_status[my_clust] = open_trx(buf,"w");
1378 clust_status_id[my_clust] = 1;
1381 else if (clust_status_id[my_clust] == -2)
1382 gmx_fatal(FARGS,"File %s.xtc should still be open (%d open xtc files)\n""in order to write frame %d. my_clust = %d",
1383 clust->grpname[my_clust],ntrxopen,frame,
1385 write_trxframe(clust_status[my_clust],&frout,gc);
1386 nfwritten[my_clust]++;
1387 if (nfwritten[my_clust] ==
1388 (clust->clust->index[my_clust+1]-
1389 clust->clust->index[my_clust])) {
1390 close_trx(clust_status[my_clust]);
1391 clust_status_id[my_clust] = -2;
1394 gmx_fatal(FARGS,"Less than zero open xtc files!");
1399 write_trxframe(trxout,&frout,gc);
1404 sprintf(title,"Generated by trjconv : %s t= %9.5f",
1406 if (bSeparate || bSplitHere)
1407 out=ffopen(out_file2,"w");
1410 write_hconf_p(out,title,&useatoms,prec2ndec(frout.prec),
1411 frout.x,fr.bV?frout.v:NULL,frout.box);
1414 fprintf(out,"REMARK GENERATED BY TRJCONV\n");
1415 sprintf(title,"%s t= %9.5f",top_title,frout.time);
1416 /* if reading from pdb, we want to keep the original
1417 model numbering else we write the output frame
1418 number plus one, because model 0 is not allowed in pdb */
1419 if (ftpin==efPDB && fr.bStep && fr.step > model_nr)
1423 write_pdbfile(out,title,&useatoms,frout.x,
1424 frout.ePBC,frout.box,' ',model_nr,gc,TRUE);
1427 frout.title = title;
1428 if (bSeparate || bTDump) {
1429 frout.bTitle = TRUE;
1431 frout.bAtoms = TRUE;
1432 frout.atoms = &useatoms;
1433 frout.bStep = FALSE;
1434 frout.bTime = FALSE;
1436 frout.bTitle = (outframe == 0);
1437 frout.bAtoms = FALSE;
1441 write_g96_conf(out,&frout,-1,NULL);
1449 gmx_fatal(FARGS,"DHE, ftp=%d\n",ftp);
1451 if (bSeparate || bSplitHere)
1454 /* execute command */
1457 sprintf(c,"%s %d",exec_command,file_nr-1);
1458 /*fprintf(stderr,"Executing '%s'\n",c);*/
1459 #ifdef GMX_NO_SYSTEM
1460 printf("Warning-- No calls to system(3) supported on this platform.");
1461 printf("Warning-- Skipping execution of 'system(\"%s\")'.", c);
1465 gmx_fatal(FARGS,"Error executing command: %s",c);
1473 bHaveNextFrame=read_next_frame(oenv,status,&fr);
1474 } while (!(bTDump && bDumpFrame) && bHaveNextFrame);
1477 if (!bHaveFirstFrame || (bTDump && !bDumpFrame))
1478 fprintf(stderr,"\nWARNING no output, "
1479 "last frame read at t=%g\n",fr.time);
1480 fprintf(stderr,"\n");
1484 gmx_rmpbc_done(gpbc);
1488 else if (out != NULL)
1491 for(i=0; (i<clust->clust->nr); i++)
1492 if (clust_status[i] )
1493 close_trx(clust_status[i]);
1497 do_view(oenv,out_file,NULL);