?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Merge branch 'master' into hbond
[gromacs/qmmm-gamess-us.git] / src / tools / gmx_trjconv.c
blob1c7a3fe41fdc4e358e8ba427c4feba2e8c4ee026
1 /*
2 *
3 * This source code is part of
4 *
5 * G R O M A C S
6 *
7 * GROningen MAchine for Chemical Simulations
8 *
9 * VERSION 3.2.0
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.
14
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.
19 *
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.
26 *
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.
29 *
30 * For more info, check our website at http://www.gromacs.org
31 *
32 * And Hey:
33 * Green Red Orange Magenta Azure Cyan Skyblue
34 */
35 #ifdef HAVE_CONFIG_H
36 #include <config.h>
37 #endif
38
39 #include <string.h>
40 #include <math.h>
41 #include "macros.h"
42 #include "assert.h"
43 #include "sysstuff.h"
44 #include "smalloc.h"
45 #include "typedefs.h"
46 #include "copyrite.h"
47 #include "gmxfio.h"
48 #include "tpxio.h"
49 #include "trnio.h"
50 #include "statutil.h"
51 #include "futil.h"
52 #include "pdbio.h"
53 #include "confio.h"
54 #include "names.h"
55 #include "index.h"
56 #include "vec.h"
57 #include "xtcio.h"
58 #include "do_fit.h"
59 #include "rmpbc.h"
60 #include "wgms.h"
61 #include "magic.h"
62 #include "pbc.h"
63 #include "viewit.h"
64 #include "xvgr.h"
65 #include "gmx_ana.h"
66
67 enum { euSel,euRect, euTric, euCompact, euNR};
68
69 static real calc_isquared(int nmol,rvec m_com[],rvec m_shift[],rvec clust_com)
70 {
71 real Isq=0;
72 int i;
73 rvec m0,dx;
74
75 for(i=0; (i<nmol); i++) {
76 rvec_add(m_com[i],m_shift[i],m0);
77 rvec_sub(m0,clust_com,dx);
78 Isq += iprod(dx,dx);
79 }
80 return Isq;
81 }
82
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)
86 {
87 const real tol=1e-3;
88 bool bChanged;
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;
93 atom_id *molind;
94 bool *bMol,*bTmp;
95 rvec *m_com,*m_shift,m0;
96 t_pbc pbc;
97
98 calc_box_center(ecenter,box,box_center);
99
100 /* Initiate the pbc structure */
101 memset(&pbc,0,sizeof(pbc));
102 set_pbc(&pbc,ePBC,box);
103
104 /* Convert atom index to molecular */
105 nmol = top->mols.nr;
106 molind = top->mols.index;
107 snew(bMol,nmol);
108 snew(m_com,nmol);
109 snew(m_shift,nmol);
110 snew(bTmp,top->atoms.nr);
111 nmol_cl = 0;
112 for(i=0; (i<nrefat); i++) {
113 /* Mark all molecules in the index */
114 ai = index[i];
115 bTmp[ai] = TRUE;
116 /* Binary search assuming the molecules are sorted */
117 j0=0;
118 j1=nmol-1;
119 while (j0 < j1) {
120 if (ai < molind[j0+1])
121 j1 = j0;
122 else if (ai >= molind[j1])
123 j0 = j1;
124 else {
125 jj = (j0+j1)/2;
126 if (ai < molind[jj+1])
127 j1 = jj;
128 else
129 j0 = jj;
130 }
131 }
132 bMol[j0] = TRUE;
133 }
134 /* Double check whether all atoms in all molecules that are marked are part
135 * of the cluster. Simultaneously compute the center of geometry.
136 */
137 min_dist2 = 10*sqr(trace(box));
138 imol_center = -1;
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",
143 i+1,j+1);
144 else if (!bMol[i] && bTmp[j])
145 gmx_fatal(FARGS,"Atom %d marked for clustering but not molecule %d",
146 j+1,i+1);
147 else if (bMol[i]) {
148 /* Compute center of geometry of molecule */
149 rvec_inc(m_com[i],x[j]);
150 }
151 }
152 if (bMol[i]) {
153 /* Normalize center of geometry */
154 fac = 1.0/(molind[i+1]-molind[i]);
155 for(m=0; (m<DIM); m++)
156 m_com[i][m] *= fac;
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);
161 imol_center = i;
162 }
163 nmol_cl++;
164 }
165 }
166 sfree(bTmp);
167
168 if (nmol_cl <= 0) {
169 fprintf(stderr,"No molecules selected in the cluster\n");
170 return;
171 } else if (imol_center == -1) {
172 fprintf(stderr,"No central molecules could be found\n");
173 return;
174 }
175
176
177 /* First calculation is incremental */
178 clear_rvec(clust_com);
179 Isq = 0;
180 for(i=m=0; (i<nmol); i++) {
181 /* Check whether this molecule is part of the cluster */
182 if (bMol[i]) {
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]);
192 }
193 else {
194 rvec_inc(clust_com,m_com[i]);
195 }
196 m++;
197 }
198 }
199 assert(m == nmol_cl);
200 svmul(1/nmol_cl,clust_com,clust_com);
201 put_atom_in_box(box,clust_com);
202
203 /* Now check if any molecule is more than half the box from the COM */
204 if (box) {
205 iter = 0;
206 do {
207 bChanged = FALSE;
208 for(i=0; (i<nmol) && !bChanged; i++) {
209 if (bMol[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;
220 bChanged = TRUE;
221 }
222 }
223 }
224 Isq = calc_isquared(nmol,m_com,m_shift,clust_com);
225 put_atom_in_box(box,clust_com);
226
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);
230 iter++;
231 } while (bChanged);
232 }
233 /* Now transfer the shift to all atoms in the molecule */
234 for(i=0; (i<nmol); i++) {
235 if (bMol[i]) {
236 for(j=molind[i]; (j<molind[i+1]); j++)
237 rvec_inc(x[j],m_shift[i]);
238 }
239 }
240 sfree(bMol);
241 sfree(m_com);
242 sfree(m_shift);
243 }
244
245 static void put_molecule_com_in_box(int unitcell_enum,int ecenter,
246 t_block *mols,
247 int natoms,t_atom atom[],
248 int ePBC,matrix box,rvec x[])
249 {
250 atom_id i,j;
251 int d;
252 rvec com,new_com,shift,dx,box_center;
253 real m;
254 double mtot;
255 t_pbc pbc;
256
257 calc_box_center(ecenter,box,box_center);
258 set_pbc(&pbc,ePBC,box);
259 if (mols->nr <= 0)
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++) {
262 /* calc COM */
263 clear_rvec(com);
264 mtot = 0;
265 for(j=mols->index[i]; (j<mols->index[i+1] && j<natoms); j++) {
266 m = atom[j].m;
267 for(d=0; d<DIM; d++)
268 com[d] += m*x[j][d];
269 mtot += m;
270 }
271 /* calculate final COM */
272 svmul(1.0/mtot, com, com);
273
274 /* check if COM is outside box */
275 copy_rvec(com,new_com);
276 switch (unitcell_enum) {
277 case euRect:
278 put_atoms_in_box(box,1,&new_com);
279 break;
280 case euTric:
281 put_atoms_in_triclinic_unitcell(ecenter,box,1,&new_com);
282 break;
283 case euCompact:
284 put_atoms_in_compact_unitcell(ePBC,ecenter,box,1,&new_com);
285 break;
286 }
287 rvec_sub(new_com,com,shift);
288 if (norm2(shift) > 0) {
289 if (debug)
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);
294 }
295 }
296 }
297 }
298
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[])
302 {
303 atom_id i, j, res_start, res_end, res_nat;
304 int d, presnr;
305 real m;
306 double mtot;
307 rvec box_center,com,new_com,shift;
308
309 calc_box_center(ecenter,box,box_center);
310
311 presnr = NOTSET;
312 res_start = 0;
313 clear_rvec(com);
314 mtot = 0;
315 for(i=0; i<natoms+1; i++) {
316 if (i == natoms || (presnr != atom[i].resind && presnr != NOTSET)) {
317 /* calculate final COM */
318 res_end = i;
319 res_nat = res_end - res_start;
320 svmul(1.0/mtot, com, com);
321
322 /* check if COM is outside box */
323 copy_rvec(com,new_com);
324 switch (unitcell_enum) {
325 case euRect:
326 put_atoms_in_box(box,1,&new_com);
327 break;
328 case euTric:
329 put_atoms_in_triclinic_unitcell(ecenter,box,1,&new_com);
330 break;
331 case euCompact:
332 put_atoms_in_compact_unitcell(ePBC,ecenter,box,1,&new_com);
333 break;
334 }
335 rvec_sub(new_com,com,shift);
336 if (norm2(shift)) {
337 if (debug)
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);
343 }
344 clear_rvec(com);
345 mtot = 0;
346
347 /* remember start of new residue */
348 res_start = i;
349 }
350 if (i < natoms) {
351 /* calc COM */
352 m = atom[i].m;
353 for(d=0; d<DIM; d++)
354 com[d] += m*x[i][d];
355 mtot += m;
356
357 presnr = atom[i].resind;
358 }
359 }
360 }
361
362 static void center_x(int ecenter,rvec x[],matrix box,int n,int nc,atom_id ci[])
363 {
364 int i,m,ai;
365 rvec cmin,cmax,box_center,dx;
366
367 if (nc > 0) {
368 copy_rvec(x[ci[0]],cmin);
369 copy_rvec(x[ci[0]],cmax);
370 for(i=0; i<nc; i++) {
371 ai=ci[i];
372 for(m=0; m<DIM; m++) {
373 if (x[ai][m] < cmin[m])
374 cmin[m]=x[ai][m];
375 else if (x[ai][m] > cmax[m])
376 cmax[m]=x[ai][m];
377 }
378 }
379 calc_box_center(ecenter,box,box_center);
380 for(m=0; m<DIM; m++)
381 dx[m] = box_center[m]-(cmin[m]+cmax[m])*0.5;
382
383 for(i=0; i<n; i++)
384 rvec_inc(x[i],dx);
385 }
386 }
387
388 static void mk_filenm(char *base,const char *ext,int ndigit,int file_nr,
389 char out_file[])
390 {
391 char nbuf[128];
392 int nd=0,fnr;
393
394 strcpy(out_file,base);
395 fnr = file_nr;
396 while(fnr > 0) {
397 fnr = fnr/10;
398 nd++;
399 }
400 if (nd < ndigit)
401 strncat(out_file,"00000000000",ndigit-nd);
402 sprintf(nbuf,"%d.",file_nr);
403 strcat(out_file,nbuf);
404 strcat(out_file,ext);
405 }
406
407 void check_trn(const char *fn)
408 {
409 if ((fn2ftp(fn) != efTRJ) && (fn2ftp(fn) != efTRR))
410 gmx_fatal(FARGS,"%s is not a trj file, exiting\n",fn);
411 }
412
413 #if (!defined WIN32 && !defined _WIN32 && !defined WIN64 && !defined _WIN64)
414 void do_trunc(const char *fn, real t0)
415 {
416 int in;
417 FILE *fp;
418 bool bStop,bOK;
419 t_trnheader sh;
420 off_t fpos;
421 char yesno[256];
422 int j;
423 real t=0;
424
425 if (t0 == -1)
426 gmx_fatal(FARGS,"You forgot to set the truncation time");
427
428 /* Check whether this is a .trj file */
429 check_trn(fn);
430
431 in = open_trn(fn,"r");
432 fp = gmx_fio_getfp(in);
433 if (fp == NULL) {
434 fprintf(stderr,"Sorry, can not trunc %s, truncation of this filetype is not supported\n",fn);
435 close_trn(in);
436 } else {
437 j = 0;
438 fpos = gmx_fio_ftell(in);
439 bStop= FALSE;
440 while (!bStop && fread_trnheader(in,&sh,&bOK)) {
441 fread_htrn(in,&sh,NULL,NULL,NULL,NULL);
442 fpos=ftell(fp);
443 t=sh.t;
444 if (t>=t0) {
445 #ifdef HAVE_FSEEKO
446 fseeko(fp,fpos,SEEK_SET);
447 #else
448 fseek(fp,fpos,SEEK_SET);
449 #endif
450 bStop=TRUE;
451 }
452 }
453 if (bStop) {
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))
458 {
459 gmx_fatal(FARGS,"Error reading user input");
460 }
461 if (strcmp(yesno,"YES") == 0) {
462 fprintf(stderr,"Once again, I'm gonna DO this...\n");
463 close_trn(in);
464 if(0 != truncate(fn,fpos))
465 {
466 gmx_fatal(FARGS,"Error truncating file %s",fn);
467 }
468 }
469 else {
470 fprintf(stderr,"Ok, I'll forget about it\n");
471 }
472 }
473 else {
474 fprintf(stderr,"Already at end of file (t=%g)...\n",t);
475 close_trn(in);
476 }
477 }
478 }
479 #endif
480
481 int gmx_trjconv(int argc,char *argv[])
482 {
483 const char *desc[] = {
484 "trjconv 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]",
502 "The program [TT]trjcat[tt] can concatenate multiple trajectory files.",
503 "[PAR]",
504 "Currently seven formats are supported for input and output:",
505 "[TT].xtc[tt], [TT].trr[tt], [TT].trj[tt], [TT].gro[tt], [TT].g96[tt],",
506 "[TT].pdb[tt] and [TT].g87[tt].",
507 "The file formats are detected from the file extension.",
508 "The precision of [TT].xtc[tt] and [TT].gro[tt] output is taken from the",
509 "input file for [TT].xtc[tt], [TT].gro[tt] and [TT].pdb[tt],",
510 "and from the [TT]-ndec[tt] option for other input formats. The precision",
511 "is always taken from [TT]-ndec[tt], when this option is set.",
512 "All other formats have fixed precision. [TT].trr[tt] and [TT].trj[tt]",
513 "output can be single or double precision, depending on the precision",
514 "of the trjconv binary.",
515 "Note that velocities are only supported in",
516 "[TT].trr[tt], [TT].trj[tt], [TT].gro[tt] and [TT].g96[tt] files.[PAR]",
517 "Option [TT]-app[tt] can be used to",
518 "append output to an existing trajectory file.",
519 "No checks are performed to ensure integrity",
520 "of the resulting combined trajectory file.[PAR]",
521 "Option [TT]-sep[tt] can be used to write every frame to a seperate",
522 ".gro, .g96 or .pdb file, default all frames all written to one file.",
523 "[TT].pdb[tt] files with all frames concatenated can be viewed with",
524 "[TT]rasmol -nmrpdb[tt].[PAR]",
525 "It is possible to select part of your trajectory and write it out",
526 "to a new trajectory file in order to save disk space, e.g. for leaving",
527 "out the water from a trajectory of a protein in water.",
528 "[BB]ALWAYS[bb] put the original trajectory on tape!",
529 "We recommend to use the portable [TT].xtc[tt] format for your analysis",
530 "to save disk space and to have portable files.[PAR]",
531 "There are two options for fitting the trajectory to a reference",
532 "either for essential dynamics analysis or for whatever.",
533 "The first option is just plain fitting to a reference structure",
534 "in the structure file, the second option is a progressive fit",
535 "in which the first timeframe is fitted to the reference structure ",
536 "in the structure file to obtain and each subsequent timeframe is ",
537 "fitted to the previously fitted structure. This way a continuous",
538 "trajectory is generated, which might not be the case when using the",
539 "regular fit method, e.g. when your protein undergoes large",
540 "conformational transitions.[PAR]",
541 "Option [TT]-pbc[tt] sets the type of periodic boundary condition",
542 "treatment:[BR]",
543 "* [TT]mol[tt] puts the center of mass of molecules in the box.[BR]",
544 "* [TT]res[tt] puts the center of mass of residues in the box.[BR]",
545 "* [TT]atom[tt] puts all the atoms in the box.[BR]",
546 "* [TT]nojump[tt] checks if atoms jump across the box and then puts",
547 "them back. This has the effect that all molecules",
548 "will remain whole (provided they were whole in the initial",
549 "conformation), note that this ensures a continuous trajectory but",
550 "molecules may diffuse out of the box. The starting configuration",
551 "for this procedure is taken from the structure file, if one is",
552 "supplied, otherwise it is the first frame.[BR]",
553 "* [TT]cluster[tt] clusters all the atoms in the selected index",
554 "such that they are all closest to the center of mass of the cluster",
555 "which is iteratively updated. Note that this will only give meaningful",
556 "results if you in fact have a cluster. Luckily that can be checked",
557 "afterwards using a trajectory viewer. Note also that if your molecules",
558 "are broken this will not work either.[BR]",
559 "* [TT]whole[tt] only makes broken molecules whole.[PAR]",
560 "Option [TT]-ur[tt] sets the unit cell representation for options",
561 "[TT]mol[tt], [TT]res[tt] and [TT]atom[tt] of [TT]-pbc[tt].",
562 "All three options give different results for triclinc boxes and",
563 "identical results for rectangular boxes.",
564 "[TT]rect[tt] is the ordinary brick shape.",
565 "[TT]tric[tt] is the triclinic unit cell.",
566 "[TT]compact[tt] puts all atoms at the closest distance from the center",
567 "of the box. This can be useful for visualizing e.g. truncated",
568 "octahedrons. The center for options [TT]tric[tt] and [TT]compact[tt]",
569 "is [TT]tric[tt] (see below), unless the option [TT]-boxcenter[tt]",
570 "is set differently.[PAR]",
571 "Option [TT]-center[tt] centers the system in the box. The user can",
572 "select the group which is used to determine the geometrical center.",
573 "Option [TT]-boxcenter[tt] sets the location of the center of the box",
574 "for options [TT]-pbc[tt] and [TT]-center[tt]. The center options are:",
575 "[TT]tric[tt]: half of the sum of the box vectors,",
576 "[TT]rect[tt]: half of the box diagonal,",
577 "[TT]zero[tt]: zero.",
578 "Use option [TT]-pbc mol[tt] in addition to [TT]-center[tt] when you",
579 "want all molecules in the box after the centering.[PAR]",
580 "With [TT]-dt[tt] it is possible to reduce the number of ",
581 "frames in the output. This option relies on the accuracy of the times",
582 "in your input trajectory, so if these are inaccurate use the",
583 "[TT]-timestep[tt] option to modify the time (this can be done",
584 "simultaneously). For making smooth movies the program [TT]g_filter[tt]",
585 "can reduce the number of frames while using low-pass frequency",
586 "filtering, this reduces aliasing of high frequency motions.[PAR]",
587 "Using [TT]-trunc[tt] trjconv can truncate [TT].trj[tt] in place, i.e.",
588 "without copying the file. This is useful when a run has crashed",
589 "during disk I/O (one more disk full), or when two contiguous",
590 "trajectories must be concatenated without have double frames.[PAR]",
591 "[TT]trjcat[tt] is more suitable for concatenating trajectory files.[PAR]",
592 "Option [TT]-dump[tt] can be used to extract a frame at or near",
593 "one specific time from your trajectory.[PAR]",
594 "Option [TT]-drop[tt] reads an [TT].xvg[tt] file with times and values.",
595 "When options [TT]-dropunder[tt] and/or [TT]-dropover[tt] are set,",
596 "frames with a value below and above the value of the respective options",
597 "will not be written."
598 };
599
600 int pbc_enum;
601 enum
602 {
603 epSel,
604 epNone,
605 epComMol,
606 epComRes,
607 epComAtom,
608 epNojump,
609 epCluster,
610 epWhole,
611 epNR
612 };
613 const char *pbc_opt[epNR + 1] =
614 { NULL, "none", "mol", "res", "atom", "nojump", "cluster", "whole",
615 NULL };
616
617 int unitcell_enum;
618 const char *unitcell_opt[euNR+1] =
619 { NULL, "rect", "tric", "compact", NULL };
620
621 enum
622 { ecSel, ecTric, ecRect, ecZero, ecNR};
623 const char *center_opt[ecNR+1] =
624 { NULL, "tric", "rect", "zero", NULL };
625 int ecenter;
626
627 int fit_enum;
628 enum
629 {
630 efSel, efNone, efFit, efFitXY, efReset, efResetXY, efPFit, efNR
631 };
632 const char *fit[efNR + 1] =
633 { NULL, "none", "rot+trans", "rotxy+transxy", "translation", "transxy",
634 "progressive", NULL };
635
636 static bool bAppend=FALSE,bSeparate=FALSE,bVels=TRUE,bForce=FALSE,bCONECT=FALSE;
637 static bool bCenter=FALSE,bTer=FALSE;
638 static int skip_nr=1,ndec=3,nzero=0;
639 static real tzero=0,delta_t=0,timestep=0,ttrunc=-1,tdump=-1,split_t=0;
640 static rvec newbox = {0,0,0}, shift = {0,0,0}, trans = {0,0,0};
641 static char *exec_command=NULL;
642 static real dropunder=0,dropover=0;
643 static bool bRound=FALSE;
644
645 t_pargs
646 pa[] =
647 {
648 { "-skip", FALSE, etINT,
649 { &skip_nr }, "Only write every nr-th frame" },
650 { "-dt", FALSE, etTIME,
651 { &delta_t },
652 "Only write frame when t MOD dt = first time (%t)" },
653 { "-round", FALSE, etBOOL,
654 { &bRound }, "Round measurements to nearest picosecond"
655 },
656 { "-dump", FALSE, etTIME,
657 { &tdump }, "Dump frame nearest specified time (%t)" },
658 { "-t0", FALSE, etTIME,
659 { &tzero },
660 "Starting time (%t) (default: don't change)" },
661 { "-timestep", FALSE, etTIME,
662 { &timestep },
663 "Change time step between input frames (%t)" },
664 { "-pbc", FALSE, etENUM,
665 { pbc_opt },
666 "PBC treatment (see help text for full description)" },
667 { "-ur", FALSE, etENUM,
668 { unitcell_opt }, "Unit-cell representation" },
669 { "-center", FALSE, etBOOL,
670 { &bCenter }, "Center atoms in box" },
671 { "-boxcenter", FALSE, etENUM,
672 { center_opt }, "Center for -pbc and -center" },
673 { "-box", FALSE, etRVEC,
674 { newbox },
675 "Size for new cubic box (default: read from input)" },
676 { "-trans", FALSE, etRVEC,
677 { trans },
678 "All coordinates will be translated by trans. This "
679 "can advantageously be combined with -pbc mol -ur "
680 "compact." },
681 { "-shift", FALSE, etRVEC,
682 { shift },
683 "All coordinates will be shifted by framenr*shift" },
684 { "-fit", FALSE, etENUM,
685 { fit },
686 "Fit molecule to ref structure in the structure file" },
687 { "-ndec", FALSE, etINT,
688 { &ndec },
689 "Precision for .xtc and .gro writing in number of "
690 "decimal places" },
691 { "-vel", FALSE, etBOOL,
692 { &bVels }, "Read and write velocities if possible" },
693 { "-force", FALSE, etBOOL,
694 { &bForce }, "Read and write forces if possible" },
695 #if (!defined WIN32 && !defined _WIN32 && !defined WIN64 && !defined _WIN64)
696 { "-trunc", FALSE, etTIME,
697 { &ttrunc },
698 "Truncate input trj file after this time (%t)" },
699 #endif
700 { "-exec", FALSE, etSTR,
701 { &exec_command },
702 "Execute command for every output frame with the "
703 "frame number as argument" },
704 { "-app", FALSE, etBOOL,
705 { &bAppend }, "Append output" },
706 { "-split", FALSE, etTIME,
707 { &split_t },
708 "Start writing new file when t MOD split = first "
709 "time (%t)" },
710 { "-sep", FALSE, etBOOL,
711 { &bSeparate },
712 "Write each frame to a separate .gro, .g96 or .pdb "
713 "file" },
714 { "-nzero", FALSE, etINT,
715 { &nzero },
716 "Prepend file number in case you use the -sep flag "
717 "with this number of zeroes" },
718 { "-ter", FALSE, etBOOL,
719 { &bTer },
720 "Use 'TER' in pdb file as end of frame in stead of "
721 "default 'ENDMDL'" },
722 { "-dropunder", FALSE, etREAL,
723 { &dropunder }, "Drop all frames below this value" },
724 { "-dropover", FALSE, etREAL,
725 { &dropover }, "Drop all frames above this value" },
726 { "-conect", FALSE, etBOOL,
727 { &bCONECT },
728 "Add conect records when writing pdb files. Useful "
729 "for visualization of non-standard molecules, e.g. "
730 "coarse grained ones" } };
731 #define NPA asize(pa)
732
733 FILE *out=NULL;
734 int trxout=NOTSET;
735 int status,ftp,ftpin=0,file_nr;
736 t_trxframe fr,frout;
737 int flags;
738 rvec *xmem=NULL,*vmem=NULL,*fmem=NULL;
739 rvec *xp=NULL,x_shift,hbox,box_center,dx;
740 real xtcpr, lambda,*w_rls=NULL;
741 int m,i,d,frame,outframe,natoms,nout,ncent,nre,newstep=0,model_nr;
742 #define SKIP 10
743 t_topology top;
744 gmx_conect gc=NULL;
745 int ePBC=-1;
746 t_atoms *atoms=NULL,useatoms;
747 matrix top_box;
748 atom_id *index,*cindex;
749 char *grpnm;
750 int *frindex,nrfri;
751 char *frname;
752 int ifit,irms,my_clust=-1;
753 atom_id *ind_fit,*ind_rms;
754 char *gn_fit,*gn_rms;
755 t_cluster_ndx *clust=NULL;
756 int *clust_status=NULL;
757 int ntrxopen=0;
758 int *nfwritten=NULL;
759 int ndrop=0,ncol,drop0=0,drop1=0,dropuse=0;
760 double **dropval;
761 real tshift=0,t0=-1,dt=0.001,prec;
762 bool bFit,bFitXY,bPFit,bReset;
763 int nfitdim;
764 bool bRmPBC,bPBCWhole,bPBCcomRes,bPBCcomMol,bPBCcomAtom,bPBC,bNoJump,bCluster;
765 bool bCopy,bDoIt,bIndex,bTDump,bSetTime,bTPS=FALSE,bDTset=FALSE;
766 bool bExec,bTimeStep=FALSE,bDumpFrame=FALSE,bSetPrec,bNeedPrec;
767 bool bHaveFirstFrame,bHaveNextFrame,bSetBox,bSetUR,bSplit=FALSE;
768 bool bSubTraj=FALSE,bDropUnder=FALSE,bDropOver=FALSE,bTrans=FALSE;
769 bool bWriteFrame,bSplitHere;
770 const char *top_file,*in_file,*out_file=NULL;
771 char out_file2[256],*charpt;
772 char *outf_base=NULL;
773 const char *outf_ext=NULL;
774 char top_title[256],title[256],command[256],filemode[5];
775 int xdr=0;
776 bool bWarnCompact=FALSE;
777 const char *warn;
778 output_env_t oenv;
779
780 t_filenm fnm[] = {
781 { efTRX, "-f", NULL, ffREAD },
782 { efTRO, "-o", NULL, ffWRITE },
783 { efTPS, NULL, NULL, ffOPTRD },
784 { efNDX, NULL, NULL, ffOPTRD },
785 { efNDX, "-fr", "frames", ffOPTRD },
786 { efNDX, "-sub", "cluster", ffOPTRD },
787 { efXVG, "-drop","drop", ffOPTRD }
788 };
789 #define NFILE asize(fnm)
790
791 CopyRight(stderr,argv[0]);
792 parse_common_args(&argc,argv,
793 PCA_CAN_BEGIN | PCA_CAN_END | PCA_CAN_VIEW |
794 PCA_TIME_UNIT | PCA_BE_NICE,
795 NFILE,fnm,NPA,pa,asize(desc),desc,
796 0,NULL,&oenv);
797
798 top_file = ftp2fn(efTPS,NFILE,fnm);
799 init_top(&top);
800
801 /* Check command line */
802 in_file=opt2fn("-f",NFILE,fnm);
803
804 if (ttrunc != -1) {
805 #if (!defined WIN32 && !defined _WIN32 && !defined WIN64 && !defined _WIN64)
806 do_trunc(in_file,ttrunc);
807 #endif
808 }
809 else {
810 /* mark active cmdline options */
811 bSetBox = opt2parg_bSet("-box", NPA, pa);
812 bSetTime = opt2parg_bSet("-t0", NPA, pa);
813 bSetPrec = opt2parg_bSet("-ndec", NPA, pa);
814 bSetUR = opt2parg_bSet("-ur", NPA, pa);
815 bExec = opt2parg_bSet("-exec", NPA, pa);
816 bTimeStep = opt2parg_bSet("-timestep", NPA, pa);
817 bTDump = opt2parg_bSet("-dump", NPA, pa);
818 bDropUnder= opt2parg_bSet("-dropunder", NPA, pa);
819 bDropOver = opt2parg_bSet("-dropover", NPA, pa);
820 bTrans = opt2parg_bSet("-trans",NPA,pa);
821 bSplit = (split_t != 0);
822
823 /* parse enum options */
824 fit_enum = nenum(fit);
825 bFit = (fit_enum==efFit || fit_enum==efFitXY);
826 bFitXY = fit_enum==efFitXY;
827 bReset = (fit_enum==efReset || fit_enum==efResetXY);
828 bPFit = fit_enum==efPFit;
829 pbc_enum = nenum(pbc_opt);
830 bPBCWhole = pbc_enum==epWhole;
831 bPBCcomRes = pbc_enum==epComRes;
832 bPBCcomMol = pbc_enum==epComMol;
833 bPBCcomAtom= pbc_enum==epComAtom;
834 bNoJump = pbc_enum==epNojump;
835 bCluster = pbc_enum==epCluster;
836 bPBC = pbc_enum!=epNone;
837 unitcell_enum = nenum(unitcell_opt);
838 ecenter = nenum(center_opt) - ecTric;
839
840 /* set and check option dependencies */
841 if (bPFit) bFit = TRUE; /* for pfit, fit *must* be set */
842 if (bFit) bReset = TRUE; /* for fit, reset *must* be set */
843 nfitdim = 0;
844 if (bFit || bReset)
845 nfitdim = (fit_enum==efFitXY || fit_enum==efResetXY) ? 2 : 3;
846 bRmPBC = bFit || bPBCWhole || bPBCcomRes || bPBCcomMol;
847 if (bSetUR) {
848 if (!(bPBCcomRes || bPBCcomMol || bPBCcomAtom)) {
849 fprintf(stderr,
850 "WARNING: Option for unitcell representation (-ur %s)\n"
851 " only has effect in combination with -pbc %s, %s or %s.\n"
852 " Ingoring unitcell representation.\n\n",
853 unitcell_opt[0],pbc_opt[2],pbc_opt[3],pbc_opt[4]);
854 bSetUR = FALSE;
855 }
856 }
857 if (bFit && bPBC) {
858 gmx_fatal(FARGS,"PBC condition treatment does not work together with rotational fit.\n"
859 "Please do the PBC condition treatment first and then run trjconv in a second step\n"
860 "for the rotational fit.\n"
861 "First doing the rotational fit and then doing the PBC treatment gives incorrect\n"
862 "results!");
863 }
864 /* set flag for pdbio to terminate frames at 'TER' (iso 'ENDMDL') */
865 pdb_use_ter(bTer);
866
867 /* ndec is in nr of decimal places, prec is a multiplication factor: */
868 prec = 1;
869 for (i=0; i<ndec; i++)
870 prec *= 10;
871
872 bIndex=ftp2bSet(efNDX,NFILE,fnm);
873
874
875 /* Determine output type */
876 out_file=opt2fn("-o",NFILE,fnm);
877 ftp=fn2ftp(out_file);
878 fprintf(stderr,"Will write %s: %s\n",ftp2ext(ftp),ftp2desc(ftp));
879 bNeedPrec = (ftp==efXTC || ftp==efGRO);
880 if (bVels) {
881 /* check if velocities are possible in input and output files */
882 ftpin=fn2ftp(in_file);
883 bVels= (ftp==efTRR || ftp==efTRJ || ftp==efGRO || ftp==efG96)
884 && (ftpin==efTRR || ftpin==efTRJ || ftpin==efGRO || ftpin==efG96 ||
885 ftpin==efCPT);
886 }
887 if (bSeparate || bSplit) {
888 outf_ext = strrchr(out_file,'.');
889 if (outf_ext == NULL)
890 gmx_fatal(FARGS,"Output file name '%s' does not contain a '.'",out_file);
891 outf_base = strdup(out_file);
892 outf_base[outf_ext - out_file] = '\0';
893 }
894
895 bSubTraj = opt2bSet("-sub",NFILE,fnm);
896 if (bSubTraj) {
897 if ((ftp != efXTC) && (ftp != efTRN))
898 gmx_fatal(FARGS,"Can only use the sub option with output file types "
899 "xtc and trr");
900 clust = cluster_index(NULL,opt2fn("-sub",NFILE,fnm));
901
902 /* Check for number of files disabled, as FOPEN_MAX is not the correct
903 * number to check for. In my linux box it is only 16.
904 */
905 if (0 && (clust->clust->nr > FOPEN_MAX-4))
906 gmx_fatal(FARGS,"Can not open enough (%d) files to write all the"
907 " trajectories.\ntry splitting the index file in %d parts.\n"
908 "FOPEN_MAX = %d",
909 clust->clust->nr,1+clust->clust->nr/FOPEN_MAX,FOPEN_MAX);
910
911 snew(clust_status,clust->clust->nr);
912 snew(nfwritten,clust->clust->nr);
913 for(i=0; (i<clust->clust->nr); i++)
914 clust_status[i] = -1;
915 bSeparate = bSplit = FALSE;
916 }
917 /* skipping */
918 if (skip_nr <= 0) {
919 }
920
921 /* Determine whether to read a topology */
922 bTPS = (ftp2bSet(efTPS,NFILE,fnm) ||
923 bRmPBC || bReset || bPBCcomMol ||
924 (ftp == efGRO) || (ftp == efPDB) || bCONECT);
925
926 /* Determine if when can read index groups */
927 bIndex = (bIndex || bTPS);
928
929 if (bTPS) {
930 read_tps_conf(top_file,top_title,&top,&ePBC,&xp,NULL,top_box,
931 bReset || bPBCcomRes);
932 atoms=&top.atoms;
933 /* top_title is only used for gro and pdb,
934 * the header in such a file is top_title t= ...
935 * to prevent a double t=, remove it from top_title
936 */
937 if ((charpt=strstr(top_title," t= ")))
938 charpt[0]='\0';
939
940 if (bCONECT)
941 gc = gmx_conect_generate(&top);
942 }
943
944 /* get frame number index */
945 frindex=NULL;
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);
949 if (debug)
950 for(i=0; i<nrfri; i++)
951 fprintf(debug,"frindex[%4d]=%4d\n",i,frindex[i]);
952 }
953
954 /* get index groups etc. */
955 if (bReset) {
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);
960
961 if (bFit) {
962 if (ifit < 2)
963 gmx_fatal(FARGS,"Need at least 2 atoms to fit!\n");
964 else if (ifit == 3)
965 fprintf(stderr,"WARNING: fitting with only 2 atoms is not unique\n");
966 }
967 }
968 else if (bCluster) {
969 printf("Select group for clustering\n");
970 get_index(atoms,ftp2fn_null(efNDX,NFILE,fnm),
971 1,&ifit,&ind_fit,&gn_fit);
972 }
973
974 if (bIndex) {
975 if (bCenter) {
976 printf("Select group for centering\n");
977 get_index(atoms,ftp2fn_null(efNDX,NFILE,fnm),
978 1,&ncent,&cindex,&grpnm);
979 }
980 printf("Select group for output\n");
981 get_index(atoms,ftp2fn_null(efNDX,NFILE,fnm),
982 1,&nout,&index,&grpnm);
983 } else {
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);
987 natoms = fr.natoms;
988 close_trj(status);
989 sfree(fr.x);
990 snew(index,natoms);
991 for(i=0;i<natoms;i++)
992 index[i]=i;
993 nout=natoms;
994 if (bCenter) {
995 ncent = nout;
996 cindex = index;
997 }
998 }
999
1000 if (bReset) {
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;
1004
1005 /* Restore reference structure and set to origin,
1006 store original location (to put structure back) */
1007 if (bRmPBC)
1008 rm_pbc(&(top.idef),ePBC,atoms->nr,top_box,xp,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]]);
1012 } else
1013 clear_rvec(x_shift);
1014
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));
1023 drop0 = 0;
1024 drop1 = 0;
1025 }
1026
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);
1037 }
1038 useatoms.nr=nout;
1039 }
1040 /* select what to read */
1041 if (ftp==efTRR || ftp==efTRJ)
1042 flags = TRX_READ_X;
1043 else
1044 flags = TRX_NEED_X;
1045 if (bVels)
1046 flags = flags | TRX_READ_V;
1047 if (bForce)
1048 flags = flags | TRX_READ_F;
1049
1050 /* open trx file for reading */
1051 bHaveFirstFrame = read_first_frame(oenv,&status,in_file,&fr,flags);
1052 if (fr.bPrec)
1053 fprintf(stderr,"\nPrecision of %s is %g (nm)\n",in_file,1/fr.prec);
1054 if (bNeedPrec) {
1055 if (bSetPrec || !fr.bPrec)
1056 fprintf(stderr,"\nSetting output precision to %g (nm)\n",1/prec);
1057 else
1058 fprintf(stderr,"Using output precision of %g (nm)\n",1/prec);
1059 }
1060
1061 if (bHaveFirstFrame) {
1062 set_trxframe_ePBC(&fr,ePBC);
1063
1064 natoms = fr.natoms;
1065
1066 if (bSetTime)
1067 tshift=tzero-fr.time;
1068 else
1069 tzero=fr.time;
1070
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);
1075 } else
1076 strcpy(filemode,"w");
1077 switch (ftp) {
1078 case efXTC:
1079 case efG87:
1080 case efTRR:
1081 case efTRJ:
1082 out=NULL;
1083 if (!bSplit && !bSubTraj)
1084 trxout = open_trx(out_file,filemode);
1085 break;
1086 case efGRO:
1087 case efG96:
1088 case efPDB:
1089 if (( !bSeparate && !bSplit ) && !bSubTraj)
1090 out=ffopen(out_file,filemode);
1091 break;
1092 }
1093
1094 bCopy = FALSE;
1095 if (bIndex)
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]);
1101 }
1102 if (bCopy) {
1103 snew(xmem,nout);
1104 if (bVels) {
1105 snew(vmem,nout);
1106 }
1107 if (bForce) {
1108 snew(fmem,nout);
1109 }
1110 }
1111
1112 if (ftp == efG87)
1113 fprintf(gmx_fio_getfp(trxout),"Generated by %s. #atoms=%d, a BOX is"
1114 " stored in this file.\n",Program(),nout);
1115
1116 /* Start the big loop over frames */
1117 file_nr = 0;
1118 frame = 0;
1119 outframe = 0;
1120 model_nr = -1;
1121 bDTset = FALSE;
1122
1123 /* Main loop over frames */
1124 do {
1125 if (!fr.bStep) {
1126 /* set the step */
1127 fr.step = newstep;
1128 newstep++;
1129 }
1130 if (bSubTraj) {
1131 /*if (frame >= clust->clust->nra)
1132 gmx_fatal(FARGS,"There are more frames in the trajectory than in the cluster index file\n");*/
1133 if (frame >= clust->maxframe)
1134 my_clust = -1;
1135 else
1136 my_clust = clust->inv_clust[frame];
1137 if ((my_clust < 0) || (my_clust >= clust->clust->nr) ||
1138 (my_clust == NO_ATID))
1139 my_clust = -1;
1140 }
1141
1142 if (bSetBox) {
1143 /* generate new box */
1144 clear_mat(fr.box);
1145 for (m=0; m<DIM; m++)
1146 fr.box[m][m] = newbox[m];
1147 }
1148
1149 if (bTrans) {
1150 for(i=0; i<natoms; i++)
1151 rvec_inc(fr.x[i],trans);
1152 }
1153
1154 if (bTDump) {
1155 /* determine timestep */
1156 if (t0 == -1) {
1157 t0 = fr.time;
1158 } else {
1159 if (!bDTset) {
1160 dt = fr.time-t0;
1161 bDTset = TRUE;
1162 }
1163 }
1164 /* This is not very elegant, as one can not dump a frame after
1165 * a timestep with is more than twice as small as the first one. */
1166 bDumpFrame = (fr.time > tdump-0.5*dt) && (fr.time <= tdump+0.5*dt);
1167 } else
1168 bDumpFrame = FALSE;
1169
1170 /* determine if an atom jumped across the box and reset it if so */
1171 if (bNoJump && (bTPS || frame!=0)) {
1172 for(d=0; d<DIM; d++)
1173 hbox[d] = 0.5*fr.box[d][d];
1174 for(i=0; i<natoms; i++) {
1175 if (bReset)
1176 rvec_dec(fr.x[i],x_shift);
1177 for(m=DIM-1; m>=0; m--)
1178 if (hbox[m] > 0) {
1179 while (fr.x[i][m]-xp[i][m] <= -hbox[m])
1180 for(d=0; d<=m; d++)
1181 fr.x[i][d] += fr.box[m][d];
1182 while (fr.x[i][m]-xp[i][m] > hbox[m])
1183 for(d=0; d<=m; d++)
1184 fr.x[i][d] -= fr.box[m][d];
1185 }
1186 }
1187 }
1188 else if (bCluster && bTPS) {
1189 rvec com;
1190
1191 calc_pbc_cluster(ecenter,ifit,&top,ePBC,fr.x,ind_fit,com,fr.box);
1192 }
1193
1194 if (bPFit) {
1195 /* Now modify the coords according to the flags,
1196 for normal fit, this is only done for output frames */
1197 if (bRmPBC)
1198 rm_pbc(&(top.idef),ePBC,natoms,fr.box,fr.x,fr.x);
1199
1200 reset_x_ndim(nfitdim,ifit,ind_fit,natoms,NULL,fr.x,w_rls);
1201 do_fit(natoms,w_rls,xp,fr.x);
1202 }
1203
1204 /* store this set of coordinates for future use */
1205 if (bPFit || bNoJump) {
1206 if (xp == NULL)
1207 snew(xp,natoms);
1208 for(i=0; (i<natoms); i++) {
1209 copy_rvec(fr.x[i],xp[i]);
1210 rvec_inc(fr.x[i],x_shift);
1211 }
1212 }
1213
1214 if ( frindex )
1215 /* see if we have a frame from the frame index group */
1216 for(i=0; i<nrfri && !bDumpFrame; i++)
1217 bDumpFrame = frame == frindex[i];
1218 if (debug && bDumpFrame)
1219 fprintf(debug,"dumping %d\n",frame);
1220
1221 bWriteFrame =
1222 ( ( !bTDump && !frindex && frame % skip_nr == 0 ) || bDumpFrame );
1223
1224 if (bWriteFrame && (bDropUnder || bDropOver)) {
1225 while (dropval[0][drop1]<fr.time && drop1+1<ndrop) {
1226 drop0 = drop1;
1227 drop1++;
1228 }
1229 if (fabs(dropval[0][drop0] - fr.time)
1230 < fabs(dropval[0][drop1] - fr.time)) {
1231 dropuse = drop0;
1232 } else {
1233 dropuse = drop1;
1234 }
1235 if ((bDropUnder && dropval[1][dropuse] < dropunder) ||
1236 (bDropOver && dropval[1][dropuse] > dropover))
1237 bWriteFrame = FALSE;
1238 }
1239
1240 if (bWriteFrame) {
1241
1242 /* calc new time */
1243 if (bTimeStep)
1244 fr.time = tzero+frame*timestep;
1245 else
1246 if (bSetTime)
1247 fr.time += tshift;
1248
1249 if (bTDump)
1250 fprintf(stderr,"\nDumping frame at t= %g %s\n",
1251 output_env_conv_time(oenv,fr.time),output_env_get_time_unit(oenv));
1252
1253 /* check for writing at each delta_t */
1254 bDoIt=(delta_t == 0);
1255 if (!bDoIt)
1256 {
1257 if (!bRound)
1258 bDoIt=bRmod(fr.time,tzero, delta_t);
1259 else
1260 /* round() is not C89 compatible, so we do this: */
1261 bDoIt=bRmod(floor(fr.time+0.5),floor(tzero+0.5),
1262 floor(delta_t+0.5));
1263 }
1264
1265 if (bDoIt || bTDump) {
1266 /* print sometimes */
1267 if ( ((outframe % SKIP) == 0) || (outframe < SKIP) )
1268 fprintf(stderr," -> frame %6d time %8.3f \r",
1269 outframe,output_env_conv_time(oenv,fr.time));
1270
1271 if (!bPFit) {
1272 /* Now modify the coords according to the flags,
1273 for PFit we did this already! */
1274
1275 if (bRmPBC)
1276 rm_pbc(&(top.idef),ePBC,natoms,fr.box,fr.x,fr.x);
1277
1278 if (bReset) {
1279 reset_x_ndim(nfitdim,ifit,ind_fit,natoms,NULL,fr.x,w_rls);
1280 if (bFit)
1281 do_fit_ndim(nfitdim,natoms,w_rls,xp,fr.x);
1282 if (!bCenter)
1283 for(i=0; i<natoms; i++)
1284 rvec_inc(fr.x[i],x_shift);
1285 }
1286
1287 if (bCenter)
1288 center_x(ecenter,fr.x,fr.box,natoms,ncent,cindex);
1289 }
1290
1291 if (bPBCcomAtom) {
1292 switch (unitcell_enum) {
1293 case euRect:
1294 put_atoms_in_box(fr.box,natoms,fr.x);
1295 break;
1296 case euTric:
1297 put_atoms_in_triclinic_unitcell(ecenter,fr.box,natoms,fr.x);
1298 break;
1299 case euCompact:
1300 warn = put_atoms_in_compact_unitcell(ePBC,ecenter,fr.box,
1301 natoms,fr.x);
1302 if (warn && !bWarnCompact) {
1303 fprintf(stderr,"\n%s\n",warn);
1304 bWarnCompact = TRUE;
1305 }
1306 break;
1307 }
1308 }
1309 if (bPBCcomRes) {
1310 put_residue_com_in_box(unitcell_enum,ecenter,
1311 natoms,atoms->atom,ePBC,fr.box,fr.x);
1312 }
1313 if (bPBCcomMol) {
1314 put_molecule_com_in_box(unitcell_enum,ecenter,
1315 &top.mols,
1316 natoms,atoms->atom,ePBC,fr.box,fr.x);
1317 }
1318 /* Copy the input trxframe struct to the output trxframe struct */
1319 frout = fr;
1320 frout.natoms = nout;
1321 if (bNeedPrec && (bSetPrec || !fr.bPrec)) {
1322 frout.bPrec = TRUE;
1323 frout.prec = prec;
1324 }
1325 if (bCopy) {
1326 frout.x = xmem;
1327 if (bVels) {
1328 frout.v = vmem;
1329 }
1330 if (bForce) {
1331 frout.f = fmem;
1332 }
1333 for(i=0; i<nout; i++) {
1334 copy_rvec(fr.x[index[i]],frout.x[i]);
1335 if (bVels && fr.bV) {
1336 copy_rvec(fr.v[index[i]],frout.v[i]);
1337 }
1338 if (bForce && fr.bF) {
1339 copy_rvec(fr.f[index[i]],frout.f[i]);
1340 }
1341 }
1342 }
1343
1344 if (opt2parg_bSet("-shift",NPA,pa))
1345 for(i=0; i<nout; i++)
1346 for (d=0; d<DIM; d++)
1347 frout.x[i][d] += outframe*shift[d];
1348
1349 if (!bRound)
1350 bSplitHere = bSplit && bRmod(fr.time,tzero, split_t);
1351 else
1352 {
1353 /* round() is not C89 compatible, so we do this: */
1354 bSplitHere = bSplit && bRmod(floor(fr.time+0.5),
1355 floor(tzero+0.5),
1356 floor(split_t+0.5));
1357 }
1358 if (bSeparate || bSplitHere)
1359 mk_filenm(outf_base,ftp2ext(ftp),nzero,file_nr,out_file2);
1360
1361 switch(ftp) {
1362 case efTRJ:
1363 case efTRR:
1364 case efG87:
1365 case efXTC:
1366 if ( bSplitHere ) {
1367 if ( trxout >= 0 )
1368 close_trx(trxout);
1369 trxout = open_trx(out_file2,filemode);
1370 }
1371 if (bSubTraj) {
1372 if (my_clust != -1) {
1373 char buf[STRLEN];
1374 if (clust_status[my_clust] == -1) {
1375 sprintf(buf,"%s.%s",clust->grpname[my_clust],ftp2ext(ftp));
1376 clust_status[my_clust] = open_trx(buf,"w");
1377 ntrxopen++;
1378 }
1379 else if (clust_status[my_clust] == -2)
1380 gmx_fatal(FARGS,"File %s.xtc should still be open (%d open xtc files)\n""in order to write frame %d. my_clust = %d",
1381 clust->grpname[my_clust],ntrxopen,frame,
1382 my_clust);
1383 write_trxframe(clust_status[my_clust],&frout,gc);
1384 nfwritten[my_clust]++;
1385 if (nfwritten[my_clust] ==
1386 (clust->clust->index[my_clust+1]-
1387 clust->clust->index[my_clust])) {
1388 close_trx(clust_status[my_clust]);
1389 clust_status[my_clust] = -2;
1390 ntrxopen--;
1391 if (ntrxopen < 0)
1392 gmx_fatal(FARGS,"Less than zero open xtc files!");
1393 }
1394 }
1395 }
1396 else
1397 write_trxframe(trxout,&frout,gc);
1398 break;
1399 case efGRO:
1400 case efG96:
1401 case efPDB:
1402 sprintf(title,"Generated by trjconv : %s t= %9.5f",
1403 top_title,fr.time);
1404 if (bSeparate || bSplitHere)
1405 out=ffopen(out_file2,"w");
1406 switch(ftp) {
1407 case efGRO:
1408 write_hconf_p(out,title,&useatoms,prec2ndec(frout.prec),
1409 frout.x,fr.bV?frout.v:NULL,frout.box);
1410 break;
1411 case efPDB:
1412 fprintf(out,"REMARK GENERATED BY TRJCONV\n");
1413 sprintf(title,"%s t= %9.5f",top_title,frout.time);
1414 /* if reading from pdb, we want to keep the original
1415 model numbering else we write the output frame
1416 number plus one, because model 0 is not allowed in pdb */
1417 if (ftpin==efPDB && fr.bStep && fr.step > model_nr)
1418 model_nr = fr.step;
1419 else
1420 model_nr++;
1421 write_pdbfile(out,title,&useatoms,frout.x,
1422 frout.ePBC,frout.box,0,model_nr,gc);
1423 break;
1424 case efG96:
1425 frout.title = title;
1426 if (bSeparate || bTDump) {
1427 frout.bTitle = TRUE;
1428 if (bTPS)
1429 frout.bAtoms = TRUE;
1430 frout.atoms = &useatoms;
1431 frout.bStep = FALSE;
1432 frout.bTime = FALSE;
1433 } else {
1434 frout.bTitle = (outframe == 0);
1435 frout.bAtoms = FALSE;
1436 frout.bStep = TRUE;
1437 frout.bTime = TRUE;
1438 }
1439 write_g96_conf(out,&frout,-1,NULL);
1440 }
1441 if (bSeparate) {
1442 ffclose(out);
1443 out = NULL;
1444 }
1445 break;
1446 default:
1447 gmx_fatal(FARGS,"DHE, ftp=%d\n",ftp);
1448 }
1449 if (bSeparate || bSplitHere)
1450 file_nr++;
1451
1452 /* execute command */
1453 if (bExec) {
1454 char c[255];
1455 sprintf(c,"%s %d",exec_command,file_nr-1);
1456 /*fprintf(stderr,"Executing '%s'\n",c);*/
1457 #ifdef GMX_NO_SYSTEM
1458 printf("Warning-- No calls to system(3) supported on this platform.");
1459 printf("Warning-- Skipping execution of 'system(\"%s\")'.", c);
1460 #else
1461 if(0 != system(c))
1462 {
1463 gmx_fatal(FARGS,"Error executing command: %s",c);
1464 }
1465 #endif
1466 }
1467 outframe++;
1468 }
1469 }
1470 frame++;
1471 bHaveNextFrame=read_next_frame(oenv,status,&fr);
1472 } while (!(bTDump && bDumpFrame) && bHaveNextFrame);
1473 }
1474
1475 if (!bHaveFirstFrame || (bTDump && !bDumpFrame))
1476 fprintf(stderr,"\nWARNING no output, "
1477 "last frame read at t=%g\n",fr.time);
1478 fprintf(stderr,"\n");
1479
1480 close_trj(status);
1481
1482 if (trxout >= 0)
1483 close_trx(trxout);
1484 else if (out != NULL)
1485 ffclose(out);
1486 if (bSubTraj) {
1487 for(i=0; (i<clust->clust->nr); i++)
1488 if (clust_status[i] >= 0)
1489 close_trx(clust_status[i]);
1490 }
1491 }
1492
1493 do_view(oenv,out_file,NULL);
1494
1495 thanx(stderr);
1496
1497 return 0;
1498 }
Atempt to write QMMM interface for Gamess-US