?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Fixed make_edi.c
[gromacs/rigid-bodies.git] / src / tools / gmx_traj.c
blob1e2e321eee373a9a2db33d0605c2399adf29c109
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 <math.h>
40 #include <string.h>
41 #include "statutil.h"
42 #include "sysstuff.h"
43 #include "typedefs.h"
44 #include "smalloc.h"
45 #include "macros.h"
46 #include "vec.h"
47 #include "pbc.h"
48 #include "copyrite.h"
49 #include "futil.h"
50 #include "statutil.h"
51 #include "index.h"
52 #include "mshift.h"
53 #include "xvgr.h"
54 #include "tpxio.h"
55 #include "rmpbc.h"
56 #include "physics.h"
57 #include "nrjac.h"
58 #include "confio.h"
59 #include "gmx_ana.h"
60
61
62 static void low_print_data(FILE *fp,real time,rvec x[],int n,atom_id *index,
63 bool bDim[])
64 {
65 int i,d;
66
67 fprintf(fp," %g",time);
68 for(i=0; i<n; i++) {
69 for(d=0; d<DIM; d++)
70 if (bDim[d]) {
71 if (index)
72 fprintf(fp,"\t%g",x[index[i]][d]);
73 else
74 fprintf(fp,"\t%g",x[i][d]);
75 }
76 if (bDim[DIM]) {
77 if (index)
78 fprintf(fp,"\t%g",norm(x[index[i]]));
79 else
80 fprintf(fp,"\t%g",norm(x[i]));
81 }
82 }
83 fprintf(fp,"\n");
84 }
85
86 static void average_data(rvec x[],rvec xav[],real *mass,
87 int ngrps,int isize[],atom_id **index)
88 {
89 int g,i,ind,d;
90 real m,mtot;
91 rvec tmp;
92 double sum[DIM];
93
94 for(g=0; g<ngrps; g++) {
95 for(d=0; d<DIM; d++)
96 sum[d] = 0;
97 clear_rvec(xav[g]);
98 mtot = 0;
99 for(i=0; i<isize[g]; i++) {
100 ind = index[g][i];
101 if (mass) {
102 m = mass[ind];
103 svmul(m,x[ind],tmp);
104 for(d=0; d<DIM; d++)
105 sum[d] += tmp[d];
106 mtot += m;
107 } else
108 for(d=0; d<DIM; d++)
109 sum[d] += x[ind][d];
110 }
111 if (mass) {
112 for(d=0; d<DIM; d++)
113 xav[g][d] = sum[d]/mtot;
114 } else {
115 /* mass=NULL, so these are forces: sum only (do not average) */
116 for(d=0; d<DIM; d++)
117 xav[g][d] = sum[d];
118 }
119 }
120 }
121
122 static void print_data(FILE *fp,real time,rvec x[],real *mass,bool bCom,
123 int ngrps,int isize[],atom_id **index,bool bDim[])
124 {
125 static rvec *xav=NULL;
126
127 if (bCom) {
128 if (xav==NULL)
129 snew(xav,ngrps);
130 average_data(x,xav,mass,ngrps,isize,index);
131 low_print_data(fp,time,xav,ngrps,NULL,bDim);
132 } else
133 low_print_data(fp,time,x,isize[0],index[0],bDim);
134 }
135
136 static void write_trx_x(int status,t_trxframe *fr,real *mass,bool bCom,
137 int ngrps,int isize[],atom_id **index)
138 {
139 static rvec *xav=NULL;
140 static t_atoms *atoms=NULL;
141 t_trxframe fr_av;
142 int i;
143
144 fr->bV = FALSE;
145 fr->bF = FALSE;
146 if (bCom) {
147 if (xav==NULL) {
148 snew(xav,ngrps);
149 snew(atoms,1);
150 *atoms = *fr->atoms;
151 snew(atoms->atom,ngrps);
152 atoms->nr = ngrps;
153 /* Note that atom and residue names will be the ones of the first atom
154 * in each group.
155 */
156 for(i=0; i<ngrps; i++) {
157 atoms->atom[i] = fr->atoms->atom[index[i][0]];
158 atoms->atomname[i] = fr->atoms->atomname[index[i][0]];
159 }
160 }
161 average_data(fr->x,xav,mass,ngrps,isize,index);
162 fr_av = *fr;
163 fr_av.natoms = ngrps;
164 fr_av.atoms = atoms;
165 fr_av.x = xav;
166 write_trxframe(status,&fr_av,NULL);
167 } else {
168 write_trxframe_indexed(status,fr,isize[0],index[0],NULL);
169 }
170 }
171
172 static void make_legend(FILE *fp,int ngrps,int isize,atom_id index[],
173 char **name,bool bCom,bool bMol,bool bDim[],
174 const output_env_t oenv)
175 {
176 char **leg;
177 const char *dimtxt[] = { " X", " Y", " Z", "" };
178 int n,i,j,d;
179
180 if (bCom)
181 n = ngrps;
182 else
183 n = isize;
184
185 snew(leg,4*n);
186 j=0;
187 for(i=0; i<n; i++)
188 for(d=0; d<=DIM; d++)
189 if (bDim[d]) {
190 snew(leg[j],STRLEN);
191 if (bMol)
192 sprintf(leg[j],"mol %d%s",index[i]+1,dimtxt[d]);
193 else if (bCom)
194 sprintf(leg[j],"%s%s",name[i],dimtxt[d]);
195 else
196 sprintf(leg[j],"atom %d%s",index[i]+1,dimtxt[d]);
197 j++;
198 }
199 xvgr_legend(fp,j,leg,oenv);
200 for(i=0; i<j; i++)
201 sfree(leg[i]);
202 sfree(leg);
203 }
204
205 static real ekrot(rvec x[],rvec v[],real mass[],int isize,atom_id index[])
206 {
207 static real **TCM=NULL,**L;
208 double tm,m0,lxx,lxy,lxz,lyy,lyz,lzz,ekrot;
209 rvec a0,ocm;
210 dvec dx,b0;
211 dvec xcm,vcm,acm;
212 int i,j,m,n;
213
214 if (TCM == NULL) {
215 snew(TCM,DIM);
216 for(i=0; i<DIM; i++)
217 snew(TCM[i],DIM);
218 snew(L,DIM);
219 for(i=0; i<DIM; i++)
220 snew(L[i],DIM);
221 }
222
223 clear_dvec(xcm);
224 clear_dvec(vcm);
225 clear_dvec(acm);
226 tm=0.0;
227 for(i=0; i<isize; i++) {
228 j = index[i];
229 m0=mass[j];
230 tm+=m0;
231 cprod(x[j],v[j],a0);
232 for(m=0; (m<DIM); m++) {
233 xcm[m]+=m0*x[j][m]; /* c.o.m. position */
234 vcm[m]+=m0*v[j][m]; /* c.o.m. velocity */
235 acm[m]+=m0*a0[m]; /* rotational velocity around c.o.m. */
236 }
237 }
238 dcprod(xcm,vcm,b0);
239 for(m=0; (m<DIM); m++) {
240 xcm[m]/=tm;
241 vcm[m]/=tm;
242 acm[m]-=b0[m]/tm;
243 }
244
245 lxx=lxy=lxz=lyy=lyz=lzz=0.0;
246 for(i=0; i<isize; i++) {
247 j = index[i];
248 m0=mass[j];
249 for(m=0; (m<DIM); m++)
250 dx[m]=x[j][m]-xcm[m];
251 lxx+=dx[XX]*dx[XX]*m0;
252 lxy+=dx[XX]*dx[YY]*m0;
253 lxz+=dx[XX]*dx[ZZ]*m0;
254 lyy+=dx[YY]*dx[YY]*m0;
255 lyz+=dx[YY]*dx[ZZ]*m0;
256 lzz+=dx[ZZ]*dx[ZZ]*m0;
257 }
258
259 L[XX][XX]=lyy+lzz;
260 L[YY][XX]=-lxy;
261 L[ZZ][XX]=-lxz;
262 L[XX][YY]=-lxy;
263 L[YY][YY]=lxx+lzz;
264 L[ZZ][YY]=-lyz;
265 L[XX][ZZ]=-lxz;
266 L[YY][ZZ]=-lyz;
267 L[ZZ][ZZ]=lxx+lyy;
268
269 m_inv_gen(L,DIM,TCM);
270
271 /* Compute omega (hoeksnelheid) */
272 clear_rvec(ocm);
273 ekrot=0;
274 for(m=0; (m<DIM); m++) {
275 for(n=0; (n<DIM); n++)
276 ocm[m]+=TCM[m][n]*acm[n];
277 ekrot+=0.5*ocm[m]*acm[m];
278 }
279 return ekrot;
280 }
281
282 static real ektrans(rvec v[],real mass[],int isize,atom_id index[])
283 {
284 dvec mvcom;
285 real mtot=0;
286 int i,j,d;
287
288 clear_dvec(mvcom);
289 for(i=0; i<isize; i++) {
290 j = index[i];
291 for(d=0; d<DIM; d++)
292 mvcom[d] += mass[j]*v[j][d];
293 mtot += mass[j];
294 }
295
296 return dnorm2(mvcom)/(mtot*2);
297 }
298
299 static real temp(rvec v[],real mass[],int isize,atom_id index[])
300 {
301 double ekin2=0;
302 int i,j;
303
304 for(i=0; i<isize; i++) {
305 j = index[i];
306 ekin2 += mass[j]*norm2(v[j]);
307 }
308
309 return ekin2/(3*isize*BOLTZ);
310 }
311
312 static void remove_jump(matrix box,int natoms,rvec xp[],rvec x[])
313 {
314 rvec hbox;
315 int d,i,m;
316
317 for(d=0; d<DIM; d++)
318 hbox[d] = 0.5*box[d][d];
319 for(i=0; i<natoms; i++)
320 for(m=DIM-1; m>=0; m--) {
321 while (x[i][m]-xp[i][m] <= -hbox[m])
322 for(d=0; d<=m; d++)
323 x[i][d] += box[m][d];
324 while (x[i][m]-xp[i][m] > hbox[m])
325 for(d=0; d<=m; d++)
326 x[i][d] -= box[m][d];
327 }
328 }
329
330 static void write_pdb_bfac(const char *fname,const char *xname,
331 const char *title,t_atoms *atoms,int ePBC,matrix box,
332 int isize,atom_id *index,int nfr_x,rvec *x,
333 int nfr_v,rvec *sum,
334 bool bDim[],real scale_factor,
335 const output_env_t oenv)
336 {
337 FILE *fp;
338 real max,len2,scale;
339 atom_id maxi;
340 int i,m,onedim;
341 bool bOne;
342
343 if ((nfr_x == 0) || (nfr_v == 0)) {
344 fprintf(stderr,"No frames found for %s, will not write %s\n",title,fname);
345 } else {
346 fprintf(stderr,"Used %d frames for %s\n",nfr_x,"coordinates");
347 fprintf(stderr,"Used %d frames for %s\n",nfr_v,title);
348 onedim = -1;
349 if (!bDim[DIM]) {
350 m = 0;
351 for(i=0; i<DIM; i++)
352 if (bDim[i]) {
353 onedim = i;
354 m++;
355 }
356 if (m != 1)
357 onedim = -1;
358 }
359 scale = 1.0/nfr_x;
360 for(i=0; i<isize; i++)
361 svmul(scale,x[index[i]],x[index[i]]);
362 scale = 1.0/nfr_v;
363 for(i=0; i<isize; i++)
364 svmul(scale,sum[index[i]],sum[index[i]]);
365
366 fp=xvgropen(xname,title,"Atom","",oenv);
367 for(i=0; i<isize; i++)
368 fprintf(fp,"%-5d %10.3f %10.3f %10.3f\n",i,
369 sum[index[i]][XX],sum[index[i]][YY],sum[index[i]][ZZ]);
370 ffclose(fp);
371 max = 0;
372 maxi = 0;
373 for(i=0; i<isize; i++) {
374 len2 = 0;
375 for(m=0; m<DIM; m++)
376 if (bDim[m] || bDim[DIM])
377 len2 += sqr(sum[index[i]][m]);
378 if (len2 > max) {
379 max = len2;
380 maxi = index[i];
381 }
382 }
383 if (scale_factor != 0) {
384 scale = scale_factor;
385 } else {
386 if (max == 0)
387 scale = 1;
388 else
389 scale = 10.0/sqrt(max);
390 }
391
392 printf("Maximum %s is %g on atom %d %s, res. %s %d\n",
393 title,sqrt(max)/nfr_v,maxi+1,*(atoms->atomname[maxi]),
394 *(atoms->resinfo[atoms->atom[maxi].resind].name),
395 atoms->resinfo[atoms->atom[maxi].resind].nr);
396
397 if (atoms->pdbinfo == NULL)
398 snew(atoms->pdbinfo,atoms->nr);
399 if (onedim == -1) {
400 for(i=0; i<isize; i++) {
401 len2 = 0;
402 for(m=0; m<DIM; m++)
403 if (bDim[m] || bDim[DIM])
404 len2 += sqr(sum[index[i]][m]);
405 atoms->pdbinfo[index[i]].bfac = sqrt(len2)*scale;
406 }
407 } else {
408 for(i=0; i<isize; i++)
409 atoms->pdbinfo[index[i]].bfac = sum[index[i]][onedim]*scale;
410 }
411 write_sto_conf_indexed(fname,title,atoms,x,NULL,ePBC,box,isize,index);
412 }
413 }
414
415 static void update_histo(int gnx,atom_id index[],rvec v[],
416 int *nhisto,int **histo,real binwidth)
417 {
418 int i,m,in,nnn;
419 real vn,vnmax;
420
421 if (histo == NULL) {
422 vnmax = 0;
423 for(i=0; (i<gnx); i++) {
424 vn = norm(v[index[i]]);
425 vnmax = max(vn,vnmax);
426 }
427 vnmax *= 2;
428 *nhisto = 1+(vnmax/binwidth);
429 snew(*histo,*nhisto);
430 }
431 for(i=0; (i<gnx); i++) {
432 vn = norm(v[index[i]]);
433 in = vn/binwidth;
434 if (in >= *nhisto) {
435 nnn = in+100;
436 fprintf(stderr,"Extending histogram from %d to %d\n",*nhisto,nnn);
437
438 srenew(*histo,nnn);
439 for(m=*nhisto; (m<nnn); m++)
440 (*histo)[m] = 0;
441 *nhisto = nnn;
442 }
443 (*histo)[in]++;
444 }
445 }
446
447 static void print_histo(const char *fn,int nhisto,int histo[],real binwidth,
448 const output_env_t oenv)
449 {
450 FILE *fp;
451 int i;
452
453 fp = xvgropen(fn,"Velocity distribution","V (nm/ps)","arbitrary units",
454 oenv);
455 for(i=0; (i<nhisto); i++)
456 fprintf(fp,"%10.3e %10d\n",i*binwidth,histo[i]);
457 ffclose(fp);
458 }
459
460 int gmx_traj(int argc,char *argv[])
461 {
462 const char *desc[] = {
463 "g_traj plots coordinates, velocities, forces and/or the box.",
464 "With [TT]-com[tt] the coordinates, velocities and forces are",
465 "calculated for the center of mass of each group.",
466 "When [TT]-mol[tt] is set, the numbers in the index file are",
467 "interpreted as molecule numbers and the same procedure as with",
468 "[TT]-com[tt] is used for each molecule.[PAR]",
469 "Option [TT]-ot[tt] plots the temperature of each group,",
470 "provided velocities are present in the trajectory file.",
471 "No corrections are made for constrained degrees of freedom!",
472 "This implies [TT]-com[tt].[PAR]",
473 "Options [TT]-ekt[tt] and [TT]-ekr[tt] plot the translational and",
474 "rotational kinetic energy of each group,",
475 "provided velocities are present in the trajectory file.",
476 "This implies [TT]-com[tt].[PAR]",
477 "Options [TT]-cv[tt] and [TT]-cf[tt] write the average velocities",
478 "and average forces as temperature factors to a pdb file with",
479 "the average coordinates. The temperature factors are scaled such",
480 "that the maximum is 10. The scaling can be changed with the option",
481 "[TT]-scale[tt]. To get the velocities or forces of one",
482 "frame set both [TT]-b[tt] and [TT]-e[tt] to the time of",
483 "desired frame. When averaging over frames you might need to use",
484 "the [TT]-nojump[tt] option to obtain the correct average coordinates.",
485 "If you select either of these option the average force and velocity",
486 "for each atom are written to an xvg file as well",
487 "(specified with [TT]-av[tt] or [TT]-af[tt]).[PAR]",
488 "Option [TT]-vd[tt] computes a velocity distribution, i.e. the",
489 "norm of the vector is plotted. In addition in the same graph",
490 "the kinetic energy distribution is given."
491 };
492 static bool bMol=FALSE,bCom=FALSE,bNoJump=FALSE;
493 static bool bX=TRUE,bY=TRUE,bZ=TRUE,bNorm=FALSE;
494 static int ngroups=1;
495 static real scale=0,binwidth=1;
496 t_pargs pa[] = {
497 { "-com", FALSE, etBOOL, {&bCom},
498 "Plot data for the com of each group" },
499 { "-mol", FALSE, etBOOL, {&bMol},
500 "Index contains molecule numbers iso atom numbers" },
501 { "-nojump", FALSE, etBOOL, {&bNoJump},
502 "Remove jumps of atoms across the box" },
503 { "-x", FALSE, etBOOL, {&bX},
504 "Plot X-component" },
505 { "-y", FALSE, etBOOL, {&bY},
506 "Plot Y-component" },
507 { "-z", FALSE, etBOOL, {&bZ},
508 "Plot Z-component" },
509 { "-ng", FALSE, etINT, {&ngroups},
510 "Number of groups to consider" },
511 { "-len", FALSE, etBOOL, {&bNorm},
512 "Plot vector length" },
513 { "-bin", FALSE, etREAL, {&binwidth},
514 "Binwidth for velocity histogram (nm/ps)" },
515 { "-scale", FALSE, etREAL, {&scale},
516 "Scale factor for pdb output, 0 is autoscale" }
517
518 };
519 FILE *outx=NULL,*outv=NULL,*outf=NULL,*outb=NULL,*outt=NULL;
520 FILE *outekt=NULL,*outekr=NULL;
521 t_topology top;
522 int ePBC;
523 real *mass,time;
524 char title[STRLEN];
525 const char *indexfn;
526 t_trxframe fr,frout;
527 int flags,nvhisto=0,*vhisto=NULL;
528 rvec *xtop,*xp=NULL;
529 rvec *sumxv=NULL,*sumv=NULL,*sumxf=NULL,*sumf=NULL;
530 matrix topbox;
531 int status,status_out=-1;
532 int i,j,n;
533 int nr_xfr,nr_vfr,nr_ffr;
534 char **grpname;
535 int *isize0,*isize;
536 atom_id **index0,**index;
537 atom_id *atndx;
538 t_block *mols;
539 bool bTop,bOX,bOXT,bOV,bOF,bOB,bOT,bEKT,bEKR,bCV,bCF;
540 bool bDim[4],bDum[4],bVD;
541 char *box_leg[6] = { "XX", "YY", "ZZ", "YX", "ZX", "ZY" };
542 output_env_t oenv;
543
544 t_filenm fnm[] = {
545 { efTRX, "-f", NULL, ffREAD },
546 { efTPS, NULL, NULL, ffREAD },
547 { efNDX, NULL, NULL, ffOPTRD },
548 { efXVG, "-ox", "coord.xvg", ffOPTWR },
549 { efTRX, "-oxt","coord.xtc", ffOPTWR },
550 { efXVG, "-ov", "veloc.xvg", ffOPTWR },
551 { efXVG, "-of", "force.xvg", ffOPTWR },
552 { efXVG, "-ob", "box.xvg", ffOPTWR },
553 { efXVG, "-ot", "temp.xvg", ffOPTWR },
554 { efXVG, "-ekt","ektrans.xvg", ffOPTWR },
555 { efXVG, "-ekr","ekrot.xvg", ffOPTWR },
556 { efXVG, "-vd", "veldist.xvg", ffOPTWR },
557 { efPDB, "-cv", "veloc.pdb", ffOPTWR },
558 { efPDB, "-cf", "force.pdb", ffOPTWR },
559 { efXVG, "-av", "all_veloc.xvg", ffOPTWR },
560 { efXVG, "-af", "all_force.xvg", ffOPTWR }
561 };
562 #define NFILE asize(fnm)
563
564 CopyRight(stderr,argv[0]);
565 parse_common_args(&argc,argv,
566 PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW | PCA_BE_NICE,
567 NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
568
569 if (bMol)
570 fprintf(stderr,"Interpreting indexfile entries as molecules.\n"
571 "Using center of mass.\n");
572
573 bOX = opt2bSet("-ox",NFILE,fnm);
574 bOXT = opt2bSet("-oxt",NFILE,fnm);
575 bOV = opt2bSet("-ov",NFILE,fnm);
576 bOF = opt2bSet("-of",NFILE,fnm);
577 bOB = opt2bSet("-ob",NFILE,fnm);
578 bOT = opt2bSet("-ot",NFILE,fnm);
579 bEKT = opt2bSet("-ekt",NFILE,fnm);
580 bEKR = opt2bSet("-ekr",NFILE,fnm);
581 bCV = opt2bSet("-cv",NFILE,fnm) || opt2bSet("-av",NFILE,fnm);
582 bCF = opt2bSet("-cf",NFILE,fnm) || opt2bSet("-af",NFILE,fnm);
583 bVD = opt2bSet("-vd",NFILE,fnm) || opt2parg_bSet("-bin",asize(pa),pa);
584 if (bMol || bOT || bEKT || bEKR)
585 bCom = TRUE;
586
587 bDim[XX] = bX;
588 bDim[YY] = bY;
589 bDim[ZZ] = bZ;
590 bDim[DIM] = bNorm;
591
592 bTop = read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,
593 &xtop,NULL,topbox,
594 bCom && (bOX || bOXT || bOV || bOT || bEKT || bEKR));
595 sfree(xtop);
596 if ((bMol || bCV || bCF) && !bTop)
597 gmx_fatal(FARGS,"Need a run input file for option -mol, -cv or -cf");
598
599 if (bMol)
600 indexfn = ftp2fn(efNDX,NFILE,fnm);
601 else
602 indexfn = ftp2fn_null(efNDX,NFILE,fnm);
603
604 if (!(bCom && !bMol))
605 ngroups = 1;
606 snew(grpname,ngroups);
607 snew(isize0,ngroups);
608 snew(index0,ngroups);
609 get_index(&(top.atoms),indexfn,ngroups,isize0,index0,grpname);
610
611 if (bMol) {
612 mols=&(top.mols);
613 atndx = mols->index;
614 ngroups = isize0[0];
615 snew(isize,ngroups);
616 snew(index,ngroups);
617 for (i=0; i<ngroups; i++) {
618 if (index0[0][i] < 0 || index0[0][i] >= mols->nr)
619 gmx_fatal(FARGS,"Molecule index (%d) is out of range (%d-%d)",
620 index0[0][i]+1,1,mols->nr);
621 isize[i] = atndx[index0[0][i]+1] - atndx[index0[0][i]];
622 snew(index[i],isize[i]);
623 for(j=0; j<isize[i]; j++)
624 index[i][j] = atndx[index0[0][i]] + j;
625 }
626 } else {
627 isize = isize0;
628 index = index0;
629 }
630 if (bCom) {
631 snew(mass,top.atoms.nr);
632 for(i=0; i<top.atoms.nr; i++)
633 mass[i] = top.atoms.atom[i].m;
634 } else
635 mass = NULL;
636
637 flags = 0;
638 if (bOX) {
639 flags = flags | TRX_READ_X;
640 outx = xvgropen(opt2fn("-ox",NFILE,fnm),
641 bCom ? "Center of mass" : "Coordinate",
642 output_env_get_xvgr_tlabel(oenv),"Coordinate (nm)",oenv);
643 make_legend(outx,ngroups,isize0[0],index0[0],grpname,bCom,bMol,bDim,oenv);
644 }
645 if (bOXT) {
646 flags = flags | TRX_READ_X;
647 status_out = open_trx(opt2fn("-oxt",NFILE,fnm),"w");
648 }
649 if (bOV) {
650 flags = flags | TRX_READ_V;
651 outv = xvgropen(opt2fn("-ov",NFILE,fnm),
652 bCom ? "Center of mass velocity" : "Velocity",
653 output_env_get_xvgr_tlabel(oenv),"Velocity (nm/ps)",oenv);
654 make_legend(outv,ngroups,isize0[0],index0[0],grpname,bCom,bMol,bDim,oenv);
655 }
656 if (bOF) {
657 flags = flags | TRX_READ_F;
658 outf = xvgropen(opt2fn("-of",NFILE,fnm),"Force",
659 output_env_get_xvgr_tlabel(oenv),"Force (kJ mol\\S-1\\N nm\\S-1\\N)",
660 oenv);
661 make_legend(outf,ngroups,isize0[0],index0[0],grpname,bCom,bMol,bDim,oenv);
662 }
663 if (bOB) {
664 outb = xvgropen(opt2fn("-ob",NFILE,fnm),"Box vector elements",
665 output_env_get_xvgr_tlabel(oenv),"(nm)",oenv);
666
667 xvgr_legend(outb,6,box_leg,oenv);
668 }
669 if (bOT) {
670 bDum[XX] = FALSE;
671 bDum[YY] = FALSE;
672 bDum[ZZ] = FALSE;
673 bDum[DIM] = TRUE;
674 flags = flags | TRX_READ_V;
675 outt = xvgropen(opt2fn("-ot",NFILE,fnm),"Temperature",
676 output_env_get_xvgr_tlabel(oenv),"(K)", oenv);
677 make_legend(outt,ngroups,isize[0],index[0],grpname,bCom,bMol,bDum,oenv);
678 }
679 if (bEKT) {
680 bDum[XX] = FALSE;
681 bDum[YY] = FALSE;
682 bDum[ZZ] = FALSE;
683 bDum[DIM] = TRUE;
684 flags = flags | TRX_READ_V;
685 outekt = xvgropen(opt2fn("-ekt",NFILE,fnm),"Center of mass translation",
686 output_env_get_xvgr_tlabel(oenv),"Energy (kJ mol\\S-1\\N)",oenv);
687 make_legend(outekt,ngroups,isize[0],index[0],grpname,bCom,bMol,bDum,oenv);
688 }
689 if (bEKR) {
690 bDum[XX] = FALSE;
691 bDum[YY] = FALSE;
692 bDum[ZZ] = FALSE;
693 bDum[DIM] = TRUE;
694 flags = flags | TRX_READ_X | TRX_READ_V;
695 outekr = xvgropen(opt2fn("-ekr",NFILE,fnm),"Center of mass rotation",
696 output_env_get_xvgr_tlabel(oenv),"Energy (kJ mol\\S-1\\N)",oenv);
697 make_legend(outekr,ngroups,isize[0],index[0],grpname,bCom,bMol,bDum,oenv);
698 }
699 if (bVD)
700 flags = flags | TRX_READ_V;
701 if (bCV)
702 flags = flags | TRX_READ_X | TRX_READ_V;
703 if (bCF)
704 flags = flags | TRX_READ_X | TRX_READ_F;
705 if ((flags == 0) && !bOB) {
706 fprintf(stderr,"Please select one or more output file options\n");
707 exit(0);
708 }
709
710 read_first_frame(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&fr,flags);
711
712 if (bCV) {
713 snew(sumxv,fr.natoms);
714 snew(sumv,fr.natoms);
715 }
716 if (bCF) {
717 snew(sumxf,fr.natoms);
718 snew(sumf,fr.natoms);
719 }
720 nr_xfr = 0;
721 nr_vfr = 0;
722 nr_ffr = 0;
723
724 do {
725 time = output_env_conv_time(oenv,fr.time);
726
727 if (fr.bX && bNoJump && fr.bBox) {
728 if (xp)
729 remove_jump(fr.box,fr.natoms,xp,fr.x);
730 else
731 snew(xp,fr.natoms);
732 for(i=0; i<fr.natoms; i++)
733 copy_rvec(fr.x[i],xp[i]);
734 }
735
736 if (fr.bX && bCom)
737 rm_pbc(&(top.idef),ePBC,fr.natoms,fr.box,fr.x,fr.x);
738
739 if (bVD && fr.bV)
740 update_histo(isize[0],index[0],fr.v,&nvhisto,&vhisto,binwidth);
741
742 if (bOX && fr.bX)
743 print_data(outx,time,fr.x,mass,bCom,ngroups,isize,index,bDim);
744 if (bOXT && fr.bX) {
745 frout = fr;
746 if (!frout.bAtoms) {
747 frout.atoms = &top.atoms;
748 frout.bAtoms = TRUE;
749 }
750 write_trx_x(status_out,&frout,mass,bCom,ngroups,isize,index);
751 }
752 if (bOV && fr.bV)
753 print_data(outv,time,fr.v,mass,bCom,ngroups,isize,index,bDim);
754 if (bOF && fr.bF)
755 print_data(outf,time,fr.f,NULL,bCom,ngroups,isize,index,bDim);
756 if (bOB && fr.bBox)
757 fprintf(outb,"\t%g\t%g\t%g\t%g\t%g\t%g\t%g\n",fr.time,
758 fr.box[XX][XX],fr.box[YY][YY],fr.box[ZZ][ZZ],
759 fr.box[YY][XX],fr.box[ZZ][XX],fr.box[ZZ][YY]);
760 if (bOT && fr.bV) {
761 fprintf(outt," %g",time);
762 for(i=0; i<ngroups; i++)
763 fprintf(outt,"\t%g",temp(fr.v,mass,isize[i],index[i]));
764 fprintf(outt,"\n");
765 }
766 if (bEKT && fr.bV) {
767 fprintf(outekt," %g",time);
768 for(i=0; i<ngroups; i++)
769 fprintf(outekt,"\t%g",ektrans(fr.v,mass,isize[i],index[i]));
770 fprintf(outekt,"\n");
771 }
772 if (bEKR && fr.bX && fr.bV) {
773 fprintf(outekr," %g",time);
774 for(i=0; i<ngroups; i++)
775 fprintf(outekr,"\t%g",ekrot(fr.x,fr.v,mass,isize[i],index[i]));
776 fprintf(outekr,"\n");
777 }
778 if (bCV) {
779 if (fr.bX) {
780 for(i=0; i<fr.natoms; i++)
781 rvec_inc(sumxv[i],fr.x[i]);
782 nr_xfr++;
783 }
784 if (fr.bV) {
785 for(i=0; i<fr.natoms; i++)
786 rvec_inc(sumv[i],fr.v[i]);
787 nr_vfr++;
788 }
789 }
790 if (bCF) {
791 if (fr.bX) {
792 for(i=0; i<fr.natoms; i++)
793 rvec_inc(sumxf[i],fr.x[i]);
794 nr_xfr++;
795 }
796 if (fr.bF) {
797 for(i=0; i<fr.natoms; i++)
798 rvec_inc(sumf[i],fr.f[i]);
799 nr_ffr++;
800 }
801 }
802
803 } while(read_next_frame(oenv,status,&fr));
804
805
806 /* clean up a bit */
807 close_trj(status);
808
809 if (bOX) ffclose(outx);
810 if (bOXT) close_trx(status_out);
811 if (bOV) ffclose(outv);
812 if (bOF) ffclose(outf);
813 if (bOB) ffclose(outb);
814 if (bOT) ffclose(outt);
815 if (bEKT) ffclose(outekt);
816 if (bEKR) ffclose(outekr);
817
818 if (bVD)
819 print_histo(opt2fn("-vd",NFILE,fnm),nvhisto,vhisto,binwidth,oenv);
820
821 if ((bCV || bCF) && (nr_vfr>1 || nr_ffr>1) && !bNoJump)
822 fprintf(stderr,"WARNING: More than one frame was used for option -cv or -cf\n"
823 "If atoms jump across the box you should use the -nojump option\n");
824
825 if (bCV)
826 write_pdb_bfac(opt2fn("-cv",NFILE,fnm),
827 opt2fn("-av",NFILE,fnm),"average velocity",&(top.atoms),
828 ePBC,topbox,isize[0],index[0],nr_xfr,sumxv,
829 nr_vfr,sumv,bDim,scale,oenv);
830 if (bCF)
831 write_pdb_bfac(opt2fn("-cf",NFILE,fnm),
832 opt2fn("-af",NFILE,fnm),"average force",&(top.atoms),
833 ePBC,topbox,isize[0],index[0],nr_xfr,sumxf,
834 nr_ffr,sumf,bDim,scale,oenv);
835
836 /* view it */
837 view_all(oenv,NFILE, fnm);
838
839 thanx(stderr);
840
841 return 0;
842 }
843