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
38 #include "gmx_header_config.h"
68 enum { euSel
,euRect
, euTric
, euCompact
, euNR
};
71 static void calc_pbc_cluster(int ecenter
,int nrefat
,t_topology
*top
,int ePBC
,
72 rvec x
[],atom_id index
[],
73 rvec clust_com
,matrix box
, rvec clustercenter
)
75 int m
,i
,j
,j0
,j1
,jj
,ai
,aj
;
78 rvec dx
,xtest
,box_center
;
90 calc_box_center(ecenter
,box
,box_center
);
92 /* Initiate the pbc structure */
93 memset(&pbc
,0,sizeof(pbc
));
94 set_pbc(&pbc
,ePBC
,box
);
96 /* Convert atom index to molecular */
98 molind
= top
->mols
.index
;
104 snew(bTmp
,top
->atoms
.nr
);
106 for(i
=0; (i
<nrefat
); i
++) {
107 /* Mark all molecules in the index */
110 /* Binary search assuming the molecules are sorted */
114 if (ai
< molind
[j0
+1])
116 else if (ai
>= molind
[j1
])
120 if (ai
< molind
[jj
+1])
128 /* Double check whether all atoms in all molecules that are marked are part
129 * of the cluster. Simultaneously compute the center of geometry.
131 min_dist2
= 10*sqr(trace(box
));
134 for(i
=0; i
<nmol
; i
++)
136 for(j
=molind
[i
]; j
<molind
[i
+1]; j
++)
138 if (bMol
[i
] && !bTmp
[j
])
140 gmx_fatal(FARGS
,"Molecule %d marked for clustering but not atom %d in it - check your index!",i
+1,j
+1);
142 else if (!bMol
[i
] && bTmp
[j
])
144 gmx_fatal(FARGS
,"Atom %d marked for clustering but not molecule %d - this is an internal error...",j
+1,i
+1);
148 /* Make molecule whole, move 2nd and higher atom to same periodicity as 1st atom in molecule */
151 pbc_dx(&pbc
,x
[j
],x
[j
-1],dx
);
152 rvec_add(x
[j
-1],dx
,x
[j
]);
154 /* Compute center of geometry of molecule - m_com[i] was zeroed when we did snew() on it! */
155 rvec_inc(m_com
[i
],x
[j
]);
160 /* Normalize center of geometry */
161 fac
= 1.0/(molind
[i
+1]-molind
[i
]);
162 for(m
=0; (m
<DIM
); m
++)
164 /* Determine which molecule is closest to the center of the box */
165 pbc_dx(&pbc
,box_center
,m_com
[i
],dx
);
168 if (tmp_r2
< min_dist2
)
173 cluster
[ncluster
++]=i
;
180 fprintf(stderr
,"No molecules selected in the cluster\n");
183 else if (imol_center
== -1)
185 fprintf(stderr
,"No central molecules could be found\n");
190 added
[nadded
++] = imol_center
;
191 bMol
[imol_center
] = FALSE
;
193 while (nadded
<ncluster
)
195 /* Find min distance between cluster molecules and those remaining to be added */
196 min_dist2
= 10*sqr(trace(box
));
199 /* Loop over added mols */
200 for(i
=0;i
<nadded
;i
++)
203 /* Loop over all mols */
204 for(j
=0;j
<ncluster
;j
++)
207 /* check those remaining to be added */
210 pbc_dx(&pbc
,m_com
[aj
],m_com
[ai
],dx
);
212 if (tmp_r2
< min_dist2
)
222 /* Add the best molecule */
223 added
[nadded
++] = jmin
;
225 /* Calculate the shift from the ai molecule */
226 pbc_dx(&pbc
,m_com
[jmin
],m_com
[imin
],dx
);
227 rvec_add(m_com
[imin
],dx
,xtest
);
228 rvec_sub(xtest
,m_com
[jmin
],m_shift
[jmin
]);
229 rvec_inc(m_com
[jmin
],m_shift
[jmin
]);
231 for(j
=molind
[jmin
]; j
<molind
[jmin
+1]; j
++)
233 rvec_inc(x
[j
],m_shift
[jmin
]);
235 fprintf(stdout
,"\rClustering iteration %d of %d...",nadded
,ncluster
);
245 fprintf(stdout
,"\n");
248 static void put_molecule_com_in_box(int unitcell_enum
,int ecenter
,
250 int natoms
,t_atom atom
[],
251 int ePBC
,matrix box
,rvec x
[])
255 rvec com
,new_com
,shift
,dx
,box_center
;
260 calc_box_center(ecenter
,box
,box_center
);
261 set_pbc(&pbc
,ePBC
,box
);
263 gmx_fatal(FARGS
,"There are no molecule descriptions. I need a .tpr file for this pbc option.");
264 for(i
=0; (i
<mols
->nr
); i
++) {
268 for(j
=mols
->index
[i
]; (j
<mols
->index
[i
+1] && j
<natoms
); j
++) {
274 /* calculate final COM */
275 svmul(1.0/mtot
, com
, com
);
277 /* check if COM is outside box */
278 copy_rvec(com
,new_com
);
279 switch (unitcell_enum
) {
281 put_atoms_in_box(ePBC
,box
,1,&new_com
);
284 put_atoms_in_triclinic_unitcell(ecenter
,box
,1,&new_com
);
287 put_atoms_in_compact_unitcell(ePBC
,ecenter
,box
,1,&new_com
);
290 rvec_sub(new_com
,com
,shift
);
291 if (norm2(shift
) > 0) {
293 fprintf (debug
,"\nShifting position of molecule %d "
294 "by %8.3f %8.3f %8.3f\n", i
+1, PR_VEC(shift
));
295 for(j
=mols
->index
[i
]; (j
<mols
->index
[i
+1] && j
<natoms
); j
++) {
296 rvec_inc(x
[j
],shift
);
302 static void put_residue_com_in_box(int unitcell_enum
,int ecenter
,
303 int natoms
,t_atom atom
[],
304 int ePBC
,matrix box
,rvec x
[])
306 atom_id i
, j
, res_start
, res_end
, res_nat
;
310 rvec box_center
,com
,new_com
,shift
;
312 calc_box_center(ecenter
,box
,box_center
);
318 for(i
=0; i
<natoms
+1; i
++) {
319 if (i
== natoms
|| (presnr
!= atom
[i
].resind
&& presnr
!= NOTSET
)) {
320 /* calculate final COM */
322 res_nat
= res_end
- res_start
;
323 svmul(1.0/mtot
, com
, com
);
325 /* check if COM is outside box */
326 copy_rvec(com
,new_com
);
327 switch (unitcell_enum
) {
329 put_atoms_in_box(ePBC
,box
,1,&new_com
);
332 put_atoms_in_triclinic_unitcell(ecenter
,box
,1,&new_com
);
335 put_atoms_in_compact_unitcell(ePBC
,ecenter
,box
,1,&new_com
);
338 rvec_sub(new_com
,com
,shift
);
341 fprintf (debug
,"\nShifting position of residue %d (atoms %u-%u) "
342 "by %g,%g,%g\n", atom
[res_start
].resind
+1,
343 res_start
+1, res_end
+1, PR_VEC(shift
));
344 for(j
=res_start
; j
<res_end
; j
++)
345 rvec_inc(x
[j
],shift
);
350 /* remember start of new residue */
360 presnr
= atom
[i
].resind
;
365 static void center_x(int ecenter
,rvec x
[],matrix box
,int n
,int nc
,atom_id ci
[])
368 rvec cmin
,cmax
,box_center
,dx
;
371 copy_rvec(x
[ci
[0]],cmin
);
372 copy_rvec(x
[ci
[0]],cmax
);
373 for(i
=0; i
<nc
; i
++) {
375 for(m
=0; m
<DIM
; m
++) {
376 if (x
[ai
][m
] < cmin
[m
])
378 else if (x
[ai
][m
] > cmax
[m
])
382 calc_box_center(ecenter
,box
,box_center
);
384 dx
[m
] = box_center
[m
]-(cmin
[m
]+cmax
[m
])*0.5;
391 static void mk_filenm(char *base
,const char *ext
,int ndigit
,int file_nr
,
397 strcpy(out_file
,base
);
405 strncat(out_file
,"00000000000",ndigit
-nd
);
406 sprintf(nbuf
,"%d.",file_nr
);
407 strcat(out_file
,nbuf
);
408 strcat(out_file
,ext
);
411 void check_trn(const char *fn
)
413 if ((fn2ftp(fn
) != efTRJ
) && (fn2ftp(fn
) != efTRR
))
414 gmx_fatal(FARGS
,"%s is not a trajectory file, exiting\n",fn
);
417 #ifndef GMX_NATIVE_WINDOWS
418 void do_trunc(const char *fn
, real t0
)
430 gmx_fatal(FARGS
,"You forgot to set the truncation time");
432 /* Check whether this is a .trj file */
435 in
= open_trn(fn
,"r");
436 fp
= gmx_fio_getfp(in
);
438 fprintf(stderr
,"Sorry, can not trunc %s, truncation of this filetype is not supported\n",fn
);
442 fpos
= gmx_fio_ftell(in
);
444 while (!bStop
&& fread_trnheader(in
,&sh
,&bOK
)) {
445 fread_htrn(in
,&sh
,NULL
,NULL
,NULL
,NULL
);
449 gmx_fseek(fp
, fpos
, SEEK_SET
);
454 fprintf(stderr
,"Do you REALLY want to truncate this trajectory (%s) at:\n"
455 "frame %d, time %g, bytes %ld ??? (type YES if so)\n",
456 fn
,j
,t
,(long int)fpos
);
457 if(1 != scanf("%s",yesno
))
459 gmx_fatal(FARGS
,"Error reading user input");
461 if (strcmp(yesno
,"YES") == 0) {
462 fprintf(stderr
,"Once again, I'm gonna DO this...\n");
464 if(0 != truncate(fn
,fpos
))
466 gmx_fatal(FARGS
,"Error truncating file %s",fn
);
470 fprintf(stderr
,"Ok, I'll forget about it\n");
474 fprintf(stderr
,"Already at end of file (t=%g)...\n",t
);
481 int gmx_trjconv(int argc
,char *argv
[])
483 const char *desc
[] = {
484 "[TT]trjconv[tt] can convert trajectory files in many ways:[BR]",
485 "[BB]1.[bb] from one format to another[BR]",
486 "[BB]2.[bb] select a subset of atoms[BR]",
487 "[BB]3.[bb] change the periodicity representation[BR]",
488 "[BB]4.[bb] keep multimeric molecules together[BR]",
489 "[BB]5.[bb] center atoms in the box[BR]",
490 "[BB]6.[bb] fit atoms to reference structure[BR]",
491 "[BB]7.[bb] reduce the number of frames[BR]",
492 "[BB]8.[bb] change the timestamps of the frames ",
493 "([TT]-t0[tt] and [TT]-timestep[tt])[BR]",
494 "[BB]9.[bb] cut the trajectory in small subtrajectories according",
495 "to information in an index file. This allows subsequent analysis of",
496 "the subtrajectories that could, for example, be the result of a",
497 "cluster analysis. Use option [TT]-sub[tt].",
498 "This assumes that the entries in the index file are frame numbers and",
499 "dumps each group in the index file to a separate trajectory file.[BR]",
500 "[BB]10.[bb] select frames within a certain range of a quantity given",
501 "in an [TT].xvg[tt] file.[PAR]",
503 "The program [TT]trjcat[tt] is better suited for concatenating multiple trajectory files.",
506 "Currently seven formats are supported for input and output:",
507 "[TT].xtc[tt], [TT].trr[tt], [TT].trj[tt], [TT].gro[tt], [TT].g96[tt],",
508 "[TT].pdb[tt] and [TT].g87[tt].",
509 "The file formats are detected from the file extension.",
510 "The precision of [TT].xtc[tt] and [TT].gro[tt] output is taken from the",
511 "input file for [TT].xtc[tt], [TT].gro[tt] and [TT].pdb[tt],",
512 "and from the [TT]-ndec[tt] option for other input formats. The precision",
513 "is always taken from [TT]-ndec[tt], when this option is set.",
514 "All other formats have fixed precision. [TT].trr[tt] and [TT].trj[tt]",
515 "output can be single or double precision, depending on the precision",
516 "of the [TT]trjconv[tt] binary.",
517 "Note that velocities are only supported in",
518 "[TT].trr[tt], [TT].trj[tt], [TT].gro[tt] and [TT].g96[tt] files.[PAR]",
520 "Option [TT]-app[tt] can be used to",
521 "append output to an existing trajectory file.",
522 "No checks are performed to ensure integrity",
523 "of the resulting combined trajectory file.[PAR]",
525 "Option [TT]-sep[tt] can be used to write every frame to a separate",
526 "[TT].gro, .g96[tt] or [TT].pdb[tt] file. By default, all frames all written to one file.",
527 "[TT].pdb[tt] files with all frames concatenated can be viewed with",
528 "[TT]rasmol -nmrpdb[tt].[PAR]",
530 "It is possible to select part of your trajectory and write it out",
531 "to a new trajectory file in order to save disk space, e.g. for leaving",
532 "out the water from a trajectory of a protein in water.",
533 "[BB]ALWAYS[bb] put the original trajectory on tape!",
534 "We recommend to use the portable [TT].xtc[tt] format for your analysis",
535 "to save disk space and to have portable files.[PAR]",
537 "There are two options for fitting the trajectory to a reference",
538 "either for essential dynamics analysis, etc.",
539 "The first option is just plain fitting to a reference structure",
540 "in the structure file. The second option is a progressive fit",
541 "in which the first timeframe is fitted to the reference structure ",
542 "in the structure file to obtain and each subsequent timeframe is ",
543 "fitted to the previously fitted structure. This way a continuous",
544 "trajectory is generated, which might not be the case when using the",
545 "regular fit method, e.g. when your protein undergoes large",
546 "conformational transitions.[PAR]",
548 "Option [TT]-pbc[tt] sets the type of periodic boundary condition",
550 "[TT]* mol[tt] puts the center of mass of molecules in the box,",
551 "and requires a run input file to be supplied with [TT]-s[tt].[BR]",
552 "[TT]* res[tt] puts the center of mass of residues in the box.[BR]",
553 "[TT]* atom[tt] puts all the atoms in the box.[BR]",
554 "[TT]* nojump[tt] checks if atoms jump across the box and then puts",
555 "them back. This has the effect that all molecules",
556 "will remain whole (provided they were whole in the initial",
557 "conformation). [BB]Note[bb] that this ensures a continuous trajectory but",
558 "molecules may diffuse out of the box. The starting configuration",
559 "for this procedure is taken from the structure file, if one is",
560 "supplied, otherwise it is the first frame.[BR]",
561 "[TT]* cluster[tt] clusters all the atoms in the selected index",
562 "such that they are all closest to the center of mass of the cluster,",
563 "which is iteratively updated. [BB]Note[bb] that this will only give meaningful",
564 "results if you in fact have a cluster. Luckily that can be checked",
565 "afterwards using a trajectory viewer. Note also that if your molecules",
566 "are broken this will not work either.[BR]",
567 "The separate option [TT]-clustercenter[tt] can be used to specify an",
568 "approximate center for the cluster. This is useful e.g. if you have",
569 "two big vesicles, and you want to maintain their relative positions.[BR]",
570 "[TT]* whole[tt] only makes broken molecules whole.[PAR]",
572 "Option [TT]-ur[tt] sets the unit cell representation for options",
573 "[TT]mol[tt], [TT]res[tt] and [TT]atom[tt] of [TT]-pbc[tt].",
574 "All three options give different results for triclinic boxes and",
575 "identical results for rectangular boxes.",
576 "[TT]rect[tt] is the ordinary brick shape.",
577 "[TT]tric[tt] is the triclinic unit cell.",
578 "[TT]compact[tt] puts all atoms at the closest distance from the center",
579 "of the box. This can be useful for visualizing e.g. truncated octahedra",
580 "or rhombic dodecahedra. The center for options [TT]tric[tt] and [TT]compact[tt]",
581 "is [TT]tric[tt] (see below), unless the option [TT]-boxcenter[tt]",
582 "is set differently.[PAR]",
584 "Option [TT]-center[tt] centers the system in the box. The user can",
585 "select the group which is used to determine the geometrical center.",
586 "Option [TT]-boxcenter[tt] sets the location of the center of the box",
587 "for options [TT]-pbc[tt] and [TT]-center[tt]. The center options are:",
588 "[TT]tric[tt]: half of the sum of the box vectors,",
589 "[TT]rect[tt]: half of the box diagonal,",
590 "[TT]zero[tt]: zero.",
591 "Use option [TT]-pbc mol[tt] in addition to [TT]-center[tt] when you",
592 "want all molecules in the box after the centering.[PAR]",
594 "It is not always possible to use combinations of [TT]-pbc[tt],",
595 "[TT]-fit[tt], [TT]-ur[tt] and [TT]-center[tt] to do exactly what",
596 "you want in one call to [TT]trjconv[tt]. Consider using multiple",
597 "calls, and check out the GROMACS website for suggestions.[PAR]",
599 "With [TT]-dt[tt], it is possible to reduce the number of ",
600 "frames in the output. This option relies on the accuracy of the times",
601 "in your input trajectory, so if these are inaccurate use the",
602 "[TT]-timestep[tt] option to modify the time (this can be done",
603 "simultaneously). For making smooth movies, the program [TT]g_filter[tt]",
604 "can reduce the number of frames while using low-pass frequency",
605 "filtering, this reduces aliasing of high frequency motions.[PAR]",
607 "Using [TT]-trunc[tt] [TT]trjconv[tt] can truncate [TT].trj[tt] in place, i.e.",
608 "without copying the file. This is useful when a run has crashed",
609 "during disk I/O (i.e. full disk), or when two contiguous",
610 "trajectories must be concatenated without having double frames.[PAR]",
612 "Option [TT]-dump[tt] can be used to extract a frame at or near",
613 "one specific time from your trajectory.[PAR]",
615 "Option [TT]-drop[tt] reads an [TT].xvg[tt] file with times and values.",
616 "When options [TT]-dropunder[tt] and/or [TT]-dropover[tt] are set,",
617 "frames with a value below and above the value of the respective options",
618 "will not be written."
634 const char *pbc_opt
[epNR
+ 1] =
635 { NULL
, "none", "mol", "res", "atom", "nojump", "cluster", "whole",
639 const char *unitcell_opt
[euNR
+1] =
640 { NULL
, "rect", "tric", "compact", NULL
};
643 { ecSel
, ecTric
, ecRect
, ecZero
, ecNR
};
644 const char *center_opt
[ecNR
+1] =
645 { NULL
, "tric", "rect", "zero", NULL
};
651 efSel
, efNone
, efFit
, efFitXY
, efReset
, efResetXY
, efPFit
, efNR
653 const char *fit
[efNR
+ 1] =
654 { NULL
, "none", "rot+trans", "rotxy+transxy", "translation", "transxy",
655 "progressive", NULL
};
657 static gmx_bool bAppend
=FALSE
,bSeparate
=FALSE
,bVels
=TRUE
,bForce
=FALSE
,bCONECT
=FALSE
;
658 static gmx_bool bCenter
=FALSE
;
659 static int skip_nr
=1,ndec
=3,nzero
=0;
660 static real tzero
=0,delta_t
=0,timestep
=0,ttrunc
=-1,tdump
=-1,split_t
=0;
661 static rvec newbox
= {0,0,0}, shift
= {0,0,0}, trans
= {0,0,0};
662 static char *exec_command
=NULL
;
663 static real dropunder
=0,dropover
=0;
664 static gmx_bool bRound
=FALSE
;
665 static rvec clustercenter
= {0,0,0};
670 { "-skip", FALSE
, etINT
,
671 { &skip_nr
}, "Only write every nr-th frame" },
672 { "-dt", FALSE
, etTIME
,
674 "Only write frame when t MOD dt = first time (%t)" },
675 { "-round", FALSE
, etBOOL
,
676 { &bRound
}, "Round measurements to nearest picosecond"
678 { "-dump", FALSE
, etTIME
,
679 { &tdump
}, "Dump frame nearest specified time (%t)" },
680 { "-t0", FALSE
, etTIME
,
682 "Starting time (%t) (default: don't change)" },
683 { "-timestep", FALSE
, etTIME
,
685 "Change time step between input frames (%t)" },
686 { "-pbc", FALSE
, etENUM
,
688 "PBC treatment (see help text for full description)" },
689 { "-ur", FALSE
, etENUM
,
690 { unitcell_opt
}, "Unit-cell representation" },
691 { "-center", FALSE
, etBOOL
,
692 { &bCenter
}, "Center atoms in box" },
693 { "-boxcenter", FALSE
, etENUM
,
694 { center_opt
}, "Center for -pbc and -center" },
695 { "-box", FALSE
, etRVEC
,
697 "Size for new cubic box (default: read from input)" },
698 { "-clustercenter", FALSE
, etRVEC
,
700 "Optional starting point for pbc cluster option" },
701 { "-trans", FALSE
, etRVEC
,
703 "All coordinates will be translated by trans. This "
704 "can advantageously be combined with -pbc mol -ur "
706 { "-shift", FALSE
, etRVEC
,
708 "All coordinates will be shifted by framenr*shift" },
709 { "-fit", FALSE
, etENUM
,
711 "Fit molecule to ref structure in the structure file" },
712 { "-ndec", FALSE
, etINT
,
714 "Precision for .xtc and .gro writing in number of "
716 { "-vel", FALSE
, etBOOL
,
717 { &bVels
}, "Read and write velocities if possible" },
718 { "-force", FALSE
, etBOOL
,
719 { &bForce
}, "Read and write forces if possible" },
720 #ifndef GMX_NATIVE_WINDOWS
721 { "-trunc", FALSE
, etTIME
,
723 "Truncate input trajectory file after this time (%t)" },
725 { "-exec", FALSE
, etSTR
,
727 "Execute command for every output frame with the "
728 "frame number as argument" },
729 { "-app", FALSE
, etBOOL
,
730 { &bAppend
}, "Append output" },
731 { "-split", FALSE
, etTIME
,
733 "Start writing new file when t MOD split = first "
735 { "-sep", FALSE
, etBOOL
,
737 "Write each frame to a separate .gro, .g96 or .pdb "
739 { "-nzero", FALSE
, etINT
,
741 "If the -sep flag is set, use these many digits "
742 "for the file numbers and prepend zeros as needed" },
743 { "-dropunder", FALSE
, etREAL
,
744 { &dropunder
}, "Drop all frames below this value" },
745 { "-dropover", FALSE
, etREAL
,
746 { &dropover
}, "Drop all frames above this value" },
747 { "-conect", FALSE
, etBOOL
,
749 "Add conect records when writing [TT].pdb[tt] files. Useful "
750 "for visualization of non-standard molecules, e.g. "
751 "coarse grained ones" } };
752 #define NPA asize(pa)
755 t_trxstatus
*trxout
=NULL
;
757 int ftp
,ftpin
=0,file_nr
;
760 rvec
*xmem
=NULL
,*vmem
=NULL
,*fmem
=NULL
;
761 rvec
*xp
=NULL
,x_shift
,hbox
,box_center
,dx
;
762 real xtcpr
, lambda
,*w_rls
=NULL
;
763 int m
,i
,d
,frame
,outframe
,natoms
,nout
,ncent
,nre
,newstep
=0,model_nr
;
768 t_atoms
*atoms
=NULL
,useatoms
;
770 atom_id
*index
,*cindex
;
774 int ifit
,irms
,my_clust
=-1;
775 atom_id
*ind_fit
,*ind_rms
;
776 char *gn_fit
,*gn_rms
;
777 t_cluster_ndx
*clust
=NULL
;
778 t_trxstatus
**clust_status
=NULL
;
779 int *clust_status_id
=NULL
;
782 int ndrop
=0,ncol
,drop0
=0,drop1
=0,dropuse
=0;
784 real tshift
=0,t0
=-1,dt
=0.001,prec
;
785 gmx_bool bFit
,bFitXY
,bPFit
,bReset
;
787 gmx_rmpbc_t gpbc
=NULL
;
788 gmx_bool bRmPBC
,bPBCWhole
,bPBCcomRes
,bPBCcomMol
,bPBCcomAtom
,bPBC
,bNoJump
,bCluster
;
789 gmx_bool bCopy
,bDoIt
,bIndex
,bTDump
,bSetTime
,bTPS
=FALSE
,bDTset
=FALSE
;
790 gmx_bool bExec
,bTimeStep
=FALSE
,bDumpFrame
=FALSE
,bSetPrec
,bNeedPrec
;
791 gmx_bool bHaveFirstFrame
,bHaveNextFrame
,bSetBox
,bSetUR
,bSplit
=FALSE
;
792 gmx_bool bSubTraj
=FALSE
,bDropUnder
=FALSE
,bDropOver
=FALSE
,bTrans
=FALSE
;
793 gmx_bool bWriteFrame
,bSplitHere
;
794 const char *top_file
,*in_file
,*out_file
=NULL
;
795 char out_file2
[256],*charpt
;
796 char *outf_base
=NULL
;
797 const char *outf_ext
=NULL
;
798 char top_title
[256],title
[256],command
[256],filemode
[5];
800 gmx_bool bWarnCompact
=FALSE
;
805 { efTRX
, "-f", NULL
, ffREAD
},
806 { efTRO
, "-o", NULL
, ffWRITE
},
807 { efTPS
, NULL
, NULL
, ffOPTRD
},
808 { efNDX
, NULL
, NULL
, ffOPTRD
},
809 { efNDX
, "-fr", "frames", ffOPTRD
},
810 { efNDX
, "-sub", "cluster", ffOPTRD
},
811 { efXVG
, "-drop","drop", ffOPTRD
}
813 #define NFILE asize(fnm)
815 CopyRight(stderr
,argv
[0]);
816 parse_common_args(&argc
,argv
,
817 PCA_CAN_BEGIN
| PCA_CAN_END
| PCA_CAN_VIEW
|
818 PCA_TIME_UNIT
| PCA_BE_NICE
,
819 NFILE
,fnm
,NPA
,pa
,asize(desc
),desc
,
822 top_file
= ftp2fn(efTPS
,NFILE
,fnm
);
825 /* Check command line */
826 in_file
=opt2fn("-f",NFILE
,fnm
);
829 #ifndef GMX_NATIVE_WINDOWS
830 do_trunc(in_file
,ttrunc
);
834 /* mark active cmdline options */
835 bSetBox
= opt2parg_bSet("-box", NPA
, pa
);
836 bSetTime
= opt2parg_bSet("-t0", NPA
, pa
);
837 bSetPrec
= opt2parg_bSet("-ndec", NPA
, pa
);
838 bSetUR
= opt2parg_bSet("-ur", NPA
, pa
);
839 bExec
= opt2parg_bSet("-exec", NPA
, pa
);
840 bTimeStep
= opt2parg_bSet("-timestep", NPA
, pa
);
841 bTDump
= opt2parg_bSet("-dump", NPA
, pa
);
842 bDropUnder
= opt2parg_bSet("-dropunder", NPA
, pa
);
843 bDropOver
= opt2parg_bSet("-dropover", NPA
, pa
);
844 bTrans
= opt2parg_bSet("-trans",NPA
,pa
);
845 bSplit
= (split_t
!= 0);
847 /* parse enum options */
848 fit_enum
= nenum(fit
);
849 bFit
= (fit_enum
==efFit
|| fit_enum
==efFitXY
);
850 bFitXY
= fit_enum
==efFitXY
;
851 bReset
= (fit_enum
==efReset
|| fit_enum
==efResetXY
);
852 bPFit
= fit_enum
==efPFit
;
853 pbc_enum
= nenum(pbc_opt
);
854 bPBCWhole
= pbc_enum
==epWhole
;
855 bPBCcomRes
= pbc_enum
==epComRes
;
856 bPBCcomMol
= pbc_enum
==epComMol
;
857 bPBCcomAtom
= pbc_enum
==epComAtom
;
858 bNoJump
= pbc_enum
==epNojump
;
859 bCluster
= pbc_enum
==epCluster
;
860 bPBC
= pbc_enum
!=epNone
;
861 unitcell_enum
= nenum(unitcell_opt
);
862 ecenter
= nenum(center_opt
) - ecTric
;
864 /* set and check option dependencies */
865 if (bPFit
) bFit
= TRUE
; /* for pfit, fit *must* be set */
866 if (bFit
) bReset
= TRUE
; /* for fit, reset *must* be set */
869 nfitdim
= (fit_enum
==efFitXY
|| fit_enum
==efResetXY
) ? 2 : 3;
870 bRmPBC
= bFit
|| bPBCWhole
|| bPBCcomRes
|| bPBCcomMol
;
873 if (!(bPBCcomRes
|| bPBCcomMol
|| bPBCcomAtom
)) {
875 "WARNING: Option for unitcell representation (-ur %s)\n"
876 " only has effect in combination with -pbc %s, %s or %s.\n"
877 " Ingoring unitcell representation.\n\n",
878 unitcell_opt
[0],pbc_opt
[2],pbc_opt
[3],pbc_opt
[4]);
883 gmx_fatal(FARGS
,"PBC condition treatment does not work together with rotational fit.\n"
884 "Please do the PBC condition treatment first and then run trjconv in a second step\n"
885 "for the rotational fit.\n"
886 "First doing the rotational fit and then doing the PBC treatment gives incorrect\n"
890 /* ndec is in nr of decimal places, prec is a multiplication factor: */
892 for (i
=0; i
<ndec
; i
++)
895 bIndex
=ftp2bSet(efNDX
,NFILE
,fnm
);
898 /* Determine output type */
899 out_file
=opt2fn("-o",NFILE
,fnm
);
900 ftp
=fn2ftp(out_file
);
901 fprintf(stderr
,"Will write %s: %s\n",ftp2ext(ftp
),ftp2desc(ftp
));
902 bNeedPrec
= (ftp
==efXTC
|| ftp
==efGRO
);
904 /* check if velocities are possible in input and output files */
905 ftpin
=fn2ftp(in_file
);
906 bVels
= (ftp
==efTRR
|| ftp
==efTRJ
|| ftp
==efGRO
|| ftp
==efG96
)
907 && (ftpin
==efTRR
|| ftpin
==efTRJ
|| ftpin
==efGRO
|| ftpin
==efG96
||
910 if (bSeparate
|| bSplit
) {
911 outf_ext
= strrchr(out_file
,'.');
912 if (outf_ext
== NULL
)
913 gmx_fatal(FARGS
,"Output file name '%s' does not contain a '.'",out_file
);
914 outf_base
= strdup(out_file
);
915 outf_base
[outf_ext
- out_file
] = '\0';
918 bSubTraj
= opt2bSet("-sub",NFILE
,fnm
);
920 if ((ftp
!= efXTC
) && (ftp
!= efTRN
))
921 gmx_fatal(FARGS
,"Can only use the sub option with output file types "
923 clust
= cluster_index(NULL
,opt2fn("-sub",NFILE
,fnm
));
925 /* Check for number of files disabled, as FOPEN_MAX is not the correct
926 * number to check for. In my linux box it is only 16.
928 if (0 && (clust
->clust
->nr
> FOPEN_MAX
-4))
929 gmx_fatal(FARGS
,"Can not open enough (%d) files to write all the"
930 " trajectories.\ntry splitting the index file in %d parts.\n"
932 clust
->clust
->nr
,1+clust
->clust
->nr
/FOPEN_MAX
,FOPEN_MAX
);
933 gmx_warning("The -sub option could require as many open output files as there are\n"
934 "index groups in the file (%d). If you get I/O errors opening new files,\n"
935 "try reducing the number of index groups in the file, and perhaps\n"
936 "using trjconv -sub several times on different chunks of your index file.\n",
939 snew(clust_status
,clust
->clust
->nr
);
940 snew(clust_status_id
,clust
->clust
->nr
);
941 snew(nfwritten
,clust
->clust
->nr
);
942 for(i
=0; (i
<clust
->clust
->nr
); i
++)
944 clust_status
[i
] = NULL
;
945 clust_status_id
[i
] = -1;
947 bSeparate
= bSplit
= FALSE
;
953 /* Determine whether to read a topology */
954 bTPS
= (ftp2bSet(efTPS
,NFILE
,fnm
) ||
955 bRmPBC
|| bReset
|| bPBCcomMol
|| bCluster
||
956 (ftp
== efGRO
) || (ftp
== efPDB
) || bCONECT
);
958 /* Determine if when can read index groups */
959 bIndex
= (bIndex
|| bTPS
);
962 read_tps_conf(top_file
,top_title
,&top
,&ePBC
,&xp
,NULL
,top_box
,
963 bReset
|| bPBCcomRes
);
966 if (0 == top
.mols
.nr
&& (bCluster
|| bPBCcomMol
))
968 gmx_fatal(FARGS
,"Option -pbc %s requires a .tpr file for the -s option",pbc_opt
[pbc_enum
]);
971 /* top_title is only used for gro and pdb,
972 * the header in such a file is top_title t= ...
973 * to prevent a double t=, remove it from top_title
975 if ((charpt
=strstr(top_title
," t= ")))
979 gc
= gmx_conect_generate(&top
);
981 gpbc
= gmx_rmpbc_init(&top
.idef
,ePBC
,top
.atoms
.nr
,top_box
);
984 /* get frame number index */
986 if (opt2bSet("-fr",NFILE
,fnm
)) {
987 printf("Select groups of frame number indices:\n");
988 rd_index(opt2fn("-fr",NFILE
,fnm
),1,&nrfri
,(atom_id
**)&frindex
,&frname
);
990 for(i
=0; i
<nrfri
; i
++)
991 fprintf(debug
,"frindex[%4d]=%4d\n",i
,frindex
[i
]);
994 /* get index groups etc. */
996 printf("Select group for %s fit\n",
997 bFit
?"least squares":"translational");
998 get_index(atoms
,ftp2fn_null(efNDX
,NFILE
,fnm
),
999 1,&ifit
,&ind_fit
,&gn_fit
);
1003 gmx_fatal(FARGS
,"Need at least 2 atoms to fit!\n");
1005 fprintf(stderr
,"WARNING: fitting with only 2 atoms is not unique\n");
1008 else if (bCluster
) {
1009 printf("Select group for clustering\n");
1010 get_index(atoms
,ftp2fn_null(efNDX
,NFILE
,fnm
),
1011 1,&ifit
,&ind_fit
,&gn_fit
);
1016 printf("Select group for centering\n");
1017 get_index(atoms
,ftp2fn_null(efNDX
,NFILE
,fnm
),
1018 1,&ncent
,&cindex
,&grpnm
);
1020 printf("Select group for output\n");
1021 get_index(atoms
,ftp2fn_null(efNDX
,NFILE
,fnm
),
1022 1,&nout
,&index
,&grpnm
);
1024 /* no index file, so read natoms from TRX */
1025 if (!read_first_frame(oenv
,&status
,in_file
,&fr
,TRX_DONT_SKIP
))
1026 gmx_fatal(FARGS
,"Could not read a frame from %s",in_file
);
1031 for(i
=0;i
<natoms
;i
++)
1041 snew(w_rls
,atoms
->nr
);
1042 for(i
=0; (i
<ifit
); i
++)
1043 w_rls
[ind_fit
[i
]]=atoms
->atom
[ind_fit
[i
]].m
;
1045 /* Restore reference structure and set to origin,
1046 store original location (to put structure back) */
1048 gmx_rmpbc(gpbc
,top
.atoms
.nr
,top_box
,xp
);
1049 copy_rvec(xp
[index
[0]],x_shift
);
1050 reset_x_ndim(nfitdim
,ifit
,ind_fit
,atoms
->nr
,NULL
,xp
,w_rls
);
1051 rvec_dec(x_shift
,xp
[index
[0]]);
1053 clear_rvec(x_shift
);
1055 if (bDropUnder
|| bDropOver
) {
1056 /* Read the .xvg file with the drop values */
1057 fprintf(stderr
,"\nReading drop file ...");
1058 ndrop
= read_xvg(opt2fn("-drop",NFILE
,fnm
),&dropval
,&ncol
);
1059 fprintf(stderr
," %d time points\n",ndrop
);
1060 if (ndrop
== 0 || ncol
< 2)
1061 gmx_fatal(FARGS
,"Found no data points in %s",
1062 opt2fn("-drop",NFILE
,fnm
));
1067 /* Make atoms struct for output in GRO or PDB files */
1068 if ((ftp
== efGRO
) || ((ftp
== efG96
) && bTPS
) || (ftp
== efPDB
)) {
1069 /* get memory for stuff to go in .pdb file */
1070 init_t_atoms(&useatoms
,atoms
->nr
,FALSE
);
1071 sfree(useatoms
.resinfo
);
1072 useatoms
.resinfo
= atoms
->resinfo
;
1073 for(i
=0;(i
<nout
);i
++) {
1074 useatoms
.atomname
[i
]=atoms
->atomname
[index
[i
]];
1075 useatoms
.atom
[i
]=atoms
->atom
[index
[i
]];
1076 useatoms
.nres
=max(useatoms
.nres
,useatoms
.atom
[i
].resind
+1);
1080 /* select what to read */
1081 if (ftp
==efTRR
|| ftp
==efTRJ
)
1086 flags
= flags
| TRX_READ_V
;
1088 flags
= flags
| TRX_READ_F
;
1090 /* open trx file for reading */
1091 bHaveFirstFrame
= read_first_frame(oenv
,&status
,in_file
,&fr
,flags
);
1093 fprintf(stderr
,"\nPrecision of %s is %g (nm)\n",in_file
,1/fr
.prec
);
1095 if (bSetPrec
|| !fr
.bPrec
)
1096 fprintf(stderr
,"\nSetting output precision to %g (nm)\n",1/prec
);
1098 fprintf(stderr
,"Using output precision of %g (nm)\n",1/prec
);
1101 if (bHaveFirstFrame
) {
1102 set_trxframe_ePBC(&fr
,ePBC
);
1107 tshift
=tzero
-fr
.time
;
1111 /* open output for writing */
1112 if ((bAppend
) && (gmx_fexist(out_file
))) {
1113 strcpy(filemode
,"a");
1114 fprintf(stderr
,"APPENDING to existing file %s\n",out_file
);
1116 strcpy(filemode
,"w");
1123 if (!bSplit
&& !bSubTraj
)
1124 trxout
= open_trx(out_file
,filemode
);
1129 if (( !bSeparate
&& !bSplit
) && !bSubTraj
)
1130 out
=ffopen(out_file
,filemode
);
1136 /* check if index is meaningful */
1137 for(i
=0; i
<nout
; i
++) {
1138 if (index
[i
] >= natoms
)
1140 "Index[%d] %d is larger than the number of atoms in the\n"
1141 "trajectory file (%d). There is a mismatch in the contents\n"
1142 "of your -f, -s and/or -n files.",i
,index
[i
]+1,natoms
);
1143 bCopy
= bCopy
|| (i
!= index
[i
]);
1156 fprintf(gmx_fio_getfp(trx_get_fileio(trxout
)),
1157 "Generated by %s. #atoms=%d, a BOX is"
1158 " stored in this file.\n",Program(),nout
);
1160 /* Start the big loop over frames */
1167 /* Main loop over frames */
1175 /*if (frame >= clust->clust->nra)
1176 gmx_fatal(FARGS,"There are more frames in the trajectory than in the cluster index file\n");*/
1177 if (frame
> clust
->maxframe
)
1180 my_clust
= clust
->inv_clust
[frame
];
1181 if ((my_clust
< 0) || (my_clust
>= clust
->clust
->nr
) ||
1182 (my_clust
== NO_ATID
))
1187 /* generate new box */
1189 for (m
=0; m
<DIM
; m
++)
1190 fr
.box
[m
][m
] = newbox
[m
];
1194 for(i
=0; i
<natoms
; i
++)
1195 rvec_inc(fr
.x
[i
],trans
);
1199 /* determine timestep */
1208 /* This is not very elegant, as one can not dump a frame after
1209 * a timestep with is more than twice as small as the first one. */
1210 bDumpFrame
= (fr
.time
> tdump
-0.5*dt
) && (fr
.time
<= tdump
+0.5*dt
);
1214 /* determine if an atom jumped across the box and reset it if so */
1215 if (bNoJump
&& (bTPS
|| frame
!=0)) {
1216 for(d
=0; d
<DIM
; d
++)
1217 hbox
[d
] = 0.5*fr
.box
[d
][d
];
1218 for(i
=0; i
<natoms
; i
++) {
1220 rvec_dec(fr
.x
[i
],x_shift
);
1221 for(m
=DIM
-1; m
>=0; m
--)
1223 while (fr
.x
[i
][m
]-xp
[i
][m
] <= -hbox
[m
])
1225 fr
.x
[i
][d
] += fr
.box
[m
][d
];
1226 while (fr
.x
[i
][m
]-xp
[i
][m
] > hbox
[m
])
1228 fr
.x
[i
][d
] -= fr
.box
[m
][d
];
1232 else if (bCluster
) {
1235 calc_pbc_cluster(ecenter
,ifit
,&top
,ePBC
,fr
.x
,ind_fit
,com
,fr
.box
,clustercenter
);
1239 /* Now modify the coords according to the flags,
1240 for normal fit, this is only done for output frames */
1242 gmx_rmpbc_trxfr(gpbc
,&fr
);
1244 reset_x_ndim(nfitdim
,ifit
,ind_fit
,natoms
,NULL
,fr
.x
,w_rls
);
1245 do_fit(natoms
,w_rls
,xp
,fr
.x
);
1248 /* store this set of coordinates for future use */
1249 if (bPFit
|| bNoJump
) {
1252 for(i
=0; (i
<natoms
); i
++) {
1253 copy_rvec(fr
.x
[i
],xp
[i
]);
1254 rvec_inc(fr
.x
[i
],x_shift
);
1259 /* see if we have a frame from the frame index group */
1260 for(i
=0; i
<nrfri
&& !bDumpFrame
; i
++)
1261 bDumpFrame
= frame
== frindex
[i
];
1262 if (debug
&& bDumpFrame
)
1263 fprintf(debug
,"dumping %d\n",frame
);
1266 ( ( !bTDump
&& !frindex
&& frame
% skip_nr
== 0 ) || bDumpFrame
);
1268 if (bWriteFrame
&& (bDropUnder
|| bDropOver
)) {
1269 while (dropval
[0][drop1
]<fr
.time
&& drop1
+1<ndrop
) {
1273 if (fabs(dropval
[0][drop0
] - fr
.time
)
1274 < fabs(dropval
[0][drop1
] - fr
.time
)) {
1279 if ((bDropUnder
&& dropval
[1][dropuse
] < dropunder
) ||
1280 (bDropOver
&& dropval
[1][dropuse
] > dropover
))
1281 bWriteFrame
= FALSE
;
1288 fr
.time
= tzero
+frame
*timestep
;
1294 fprintf(stderr
,"\nDumping frame at t= %g %s\n",
1295 output_env_conv_time(oenv
,fr
.time
),output_env_get_time_unit(oenv
));
1297 /* check for writing at each delta_t */
1298 bDoIt
=(delta_t
== 0);
1302 bDoIt
=bRmod(fr
.time
,tzero
, delta_t
);
1304 /* round() is not C89 compatible, so we do this: */
1305 bDoIt
=bRmod(floor(fr
.time
+0.5),floor(tzero
+0.5),
1306 floor(delta_t
+0.5));
1309 if (bDoIt
|| bTDump
) {
1310 /* print sometimes */
1311 if ( ((outframe
% SKIP
) == 0) || (outframe
< SKIP
) )
1312 fprintf(stderr
," -> frame %6d time %8.3f \r",
1313 outframe
,output_env_conv_time(oenv
,fr
.time
));
1316 /* Now modify the coords according to the flags,
1317 for PFit we did this already! */
1320 gmx_rmpbc_trxfr(gpbc
,&fr
);
1323 reset_x_ndim(nfitdim
,ifit
,ind_fit
,natoms
,NULL
,fr
.x
,w_rls
);
1325 do_fit_ndim(nfitdim
,natoms
,w_rls
,xp
,fr
.x
);
1327 for(i
=0; i
<natoms
; i
++)
1328 rvec_inc(fr
.x
[i
],x_shift
);
1332 center_x(ecenter
,fr
.x
,fr
.box
,natoms
,ncent
,cindex
);
1336 switch (unitcell_enum
) {
1338 put_atoms_in_box(ePBC
,fr
.box
,natoms
,fr
.x
);
1341 put_atoms_in_triclinic_unitcell(ecenter
,fr
.box
,natoms
,fr
.x
);
1344 warn
= put_atoms_in_compact_unitcell(ePBC
,ecenter
,fr
.box
,
1346 if (warn
&& !bWarnCompact
) {
1347 fprintf(stderr
,"\n%s\n",warn
);
1348 bWarnCompact
= TRUE
;
1354 put_residue_com_in_box(unitcell_enum
,ecenter
,
1355 natoms
,atoms
->atom
,ePBC
,fr
.box
,fr
.x
);
1358 put_molecule_com_in_box(unitcell_enum
,ecenter
,
1360 natoms
,atoms
->atom
,ePBC
,fr
.box
,fr
.x
);
1362 /* Copy the input trxframe struct to the output trxframe struct */
1364 frout
.bV
= (frout
.bV
&& bVels
);
1365 frout
.bF
= (frout
.bF
&& bForce
);
1366 frout
.natoms
= nout
;
1367 if (bNeedPrec
&& (bSetPrec
|| !fr
.bPrec
)) {
1379 for(i
=0; i
<nout
; i
++) {
1380 copy_rvec(fr
.x
[index
[i
]],frout
.x
[i
]);
1382 copy_rvec(fr
.v
[index
[i
]],frout
.v
[i
]);
1385 copy_rvec(fr
.f
[index
[i
]],frout
.f
[i
]);
1390 if (opt2parg_bSet("-shift",NPA
,pa
))
1391 for(i
=0; i
<nout
; i
++)
1392 for (d
=0; d
<DIM
; d
++)
1393 frout
.x
[i
][d
] += outframe
*shift
[d
];
1396 bSplitHere
= bSplit
&& bRmod(fr
.time
,tzero
, split_t
);
1399 /* round() is not C89 compatible, so we do this: */
1400 bSplitHere
= bSplit
&& bRmod(floor(fr
.time
+0.5),
1402 floor(split_t
+0.5));
1404 if (bSeparate
|| bSplitHere
)
1405 mk_filenm(outf_base
,ftp2ext(ftp
),nzero
,file_nr
,out_file2
);
1415 trxout
= open_trx(out_file2
,filemode
);
1418 if (my_clust
!= -1) {
1420 if (clust_status_id
[my_clust
] == -1) {
1421 sprintf(buf
,"%s.%s",clust
->grpname
[my_clust
],ftp2ext(ftp
));
1422 clust_status
[my_clust
] = open_trx(buf
,"w");
1423 clust_status_id
[my_clust
] = 1;
1426 else if (clust_status_id
[my_clust
] == -2)
1427 gmx_fatal(FARGS
,"File %s.xtc should still be open (%d open .xtc files)\n""in order to write frame %d. my_clust = %d",
1428 clust
->grpname
[my_clust
],ntrxopen
,frame
,
1430 write_trxframe(clust_status
[my_clust
],&frout
,gc
);
1431 nfwritten
[my_clust
]++;
1432 if (nfwritten
[my_clust
] ==
1433 (clust
->clust
->index
[my_clust
+1]-
1434 clust
->clust
->index
[my_clust
])) {
1435 close_trx(clust_status
[my_clust
]);
1436 clust_status
[my_clust
] = NULL
;
1437 clust_status_id
[my_clust
] = -2;
1440 gmx_fatal(FARGS
,"Less than zero open .xtc files!");
1445 write_trxframe(trxout
,&frout
,gc
);
1450 sprintf(title
,"Generated by trjconv : %s t= %9.5f",
1452 if (bSeparate
|| bSplitHere
)
1453 out
=ffopen(out_file2
,"w");
1456 write_hconf_p(out
,title
,&useatoms
,prec2ndec(frout
.prec
),
1457 frout
.x
,frout
.bV
?frout
.v
:NULL
,frout
.box
);
1460 fprintf(out
,"REMARK GENERATED BY TRJCONV\n");
1461 sprintf(title
,"%s t= %9.5f",top_title
,frout
.time
);
1462 /* if reading from pdb, we want to keep the original
1463 model numbering else we write the output frame
1464 number plus one, because model 0 is not allowed in pdb */
1465 if (ftpin
==efPDB
&& fr
.bStep
&& fr
.step
> model_nr
)
1469 write_pdbfile(out
,title
,&useatoms
,frout
.x
,
1470 frout
.ePBC
,frout
.box
,' ',model_nr
,gc
,TRUE
);
1473 frout
.title
= title
;
1474 if (bSeparate
|| bTDump
) {
1475 frout
.bTitle
= TRUE
;
1477 frout
.bAtoms
= TRUE
;
1478 frout
.atoms
= &useatoms
;
1479 frout
.bStep
= FALSE
;
1480 frout
.bTime
= FALSE
;
1482 frout
.bTitle
= (outframe
== 0);
1483 frout
.bAtoms
= FALSE
;
1487 write_g96_conf(out
,&frout
,-1,NULL
);
1495 gmx_fatal(FARGS
,"DHE, ftp=%d\n",ftp
);
1497 if (bSeparate
|| bSplitHere
)
1500 /* execute command */
1503 sprintf(c
,"%s %d",exec_command
,file_nr
-1);
1504 /*fprintf(stderr,"Executing '%s'\n",c);*/
1505 #ifdef GMX_NO_SYSTEM
1506 printf("Warning-- No calls to system(3) supported on this platform.");
1507 printf("Warning-- Skipping execution of 'system(\"%s\")'.", c
);
1511 gmx_fatal(FARGS
,"Error executing command: %s",c
);
1519 bHaveNextFrame
=read_next_frame(oenv
,status
,&fr
);
1520 } while (!(bTDump
&& bDumpFrame
) && bHaveNextFrame
);
1523 if (!bHaveFirstFrame
|| (bTDump
&& !bDumpFrame
))
1524 fprintf(stderr
,"\nWARNING no output, "
1525 "last frame read at t=%g\n",fr
.time
);
1526 fprintf(stderr
,"\n");
1532 gmx_rmpbc_done(gpbc
);
1536 else if (out
!= NULL
)
1539 for(i
=0; (i
<clust
->clust
->nr
); i
++)
1540 if (clust_status
[i
] )
1541 close_trx(clust_status
[i
]);
1545 do_view(oenv
,out_file
,NULL
);