?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Partial commit of the project to remove all static variables.
[gromacs.git] / src / tools / g_anaeig.c
blobfc27b251cc5a7c85476499012a1fd2ea7e643c80
1 /*
2 * $Id$
3 *
4 * This source code is part of
5 *
6 * G R O M A C S
7 *
8 * GROningen MAchine for Chemical Simulations
9 *
10 * VERSION 3.1
11 * Copyright (c) 1991-2001, University of Groningen, The Netherlands
12 * This program is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU General Public License
14 * as published by the Free Software Foundation; either version 2
15 * of the License, or (at your option) any later version.
16 *
17 * If you want to redistribute modifications, please consider that
18 * scientific software is very special. Version control is crucial -
19 * bugs must be traceable. We will be happy to consider code for
20 * inclusion in the official distribution, but derived work must not
21 * be called official GROMACS. Details are found in the README & COPYING
22 * files - if they are missing, get the official version at www.gromacs.org.
23 *
24 * To help us fund GROMACS development, we humbly ask that you cite
25 * the papers on the package - you can find them in the top README file.
26 *
27 * For more info, check our website at http://www.gromacs.org
28 *
29 * And Hey:
30 * Gromacs Runs One Microsecond At Cannonball Speeds
31 */
32
33 #include <math.h>
34 #include <string.h>
35 #include "statutil.h"
36 #include "sysstuff.h"
37 #include "typedefs.h"
38 #include "smalloc.h"
39 #include "macros.h"
40 #include "fatal.h"
41 #include "vec.h"
42 #include "pbc.h"
43 #include "copyrite.h"
44 #include "futil.h"
45 #include "statutil.h"
46 #include "rdgroup.h"
47 #include "pdbio.h"
48 #include "confio.h"
49 #include "tpxio.h"
50 #include "matio.h"
51 #include "mshift.h"
52 #include "xvgr.h"
53 #include "do_fit.h"
54 #include "rmpbc.h"
55 #include "txtdump.h"
56 #include "eigio.h"
57
58 char *proj_unit;
59
60 static real tick_spacing(real range,int minticks)
61 {
62 real sp;
63
64 sp = 0.2*exp(log(10)*ceil(log(range)/log(10)));
65 while (range/sp < minticks-1)
66 sp = sp/2;
67
68 return sp;
69 }
70
71 static void write_xvgr_graphs(char *file,int ngraphs,
72 char *title,char *xlabel,char **ylabel,
73 int n,real *x, real **y,
74 real scale_x,bool bZero, bool bSplit)
75 {
76 FILE *out;
77 int g,i;
78 real min,max,xsp,ysp;
79
80 out=ffopen(file,"w");
81 for(g=0; g<ngraphs; g++) {
82 min=y[g][0];
83 max=y[g][0];
84 for(i=0; i<n; i++) {
85 if (y[g][i]<min) min=y[g][i];
86 if (y[g][i]>max) max=y[g][i];
87 }
88 if (bZero)
89 min=0;
90 else
91 min=min-0.1*(max-min);
92 max=max+0.1*(max-min);
93 xsp=tick_spacing((x[n-1]-x[0])*scale_x,4);
94 ysp=tick_spacing(max-min,3);
95 fprintf(out,"@ with g%d\n@ g%d on\n",g,g);
96 fprintf(out,"@ g%d autoscale type AUTO\n",g);
97 if (g==0)
98 fprintf(out,"@ title \"%s\"\n",title);
99 if (g==ngraphs-1)
100 fprintf(out,"@ xaxis label \"%s\"\n",xlabel);
101 else
102 fprintf(out,"@ xaxis ticklabel off\n");
103 fprintf(out,"@ world xmin %g\n",x[0]*scale_x);
104 fprintf(out,"@ world xmax %g\n",x[n-1]*scale_x);
105 fprintf(out,"@ world ymin %g\n",min);
106 fprintf(out,"@ world ymax %g\n",max);
107 fprintf(out,"@ view xmin 0.15\n");
108 fprintf(out,"@ view xmax 0.85\n");
109 fprintf(out,"@ view ymin %g\n",0.15+(ngraphs-1-g)*0.7/ngraphs);
110 fprintf(out,"@ view ymax %g\n",0.15+(ngraphs-g)*0.7/ngraphs);
111 fprintf(out,"@ yaxis label \"%s\"\n",ylabel[g]);
112 fprintf(out,"@ xaxis tick major %g\n",xsp);
113 fprintf(out,"@ xaxis tick minor %g\n",xsp/2);
114 fprintf(out,"@ xaxis ticklabel start type spec\n");
115 fprintf(out,"@ xaxis ticklabel start %g\n",ceil(min/xsp)*xsp);
116 fprintf(out,"@ yaxis tick major %g\n",ysp);
117 fprintf(out,"@ yaxis tick minor %g\n",ysp/2);
118 fprintf(out,"@ yaxis ticklabel start type spec\n");
119 fprintf(out,"@ yaxis ticklabel start %g\n",ceil(min/ysp)*ysp);
120 if ((min<0) && (max>0)) {
121 fprintf(out,"@ zeroxaxis bar on\n");
122 fprintf(out,"@ zeroxaxis bar linestyle 3\n");
123 }
124 for(i=0; i<n; i++) {
125 if ( bSplit && i>0 && abs(x[i])<1e-5 )
126 fprintf(out,"&\n");
127 fprintf(out,"%10.4f %10.5f\n",x[i]*scale_x,y[g][i]);
128 }
129 fprintf(out,"&\n");
130 }
131 fclose(out);
132 }
133
134 static void compare(int natoms,int n1,rvec **eigvec1,int n2,rvec **eigvec2,
135 char *Eig1File,char *Eig2File)
136 {
137 int n,neig1,neig2,nrow;
138 real **eig1,**eig2;
139 int i,j,k;
140 double sum1,sum2,trace1,trace2,sab,samsb2,tmp,ip;
141
142 n = min(n1,n2);
143
144 neig1 = read_xvg(Eig1File,&eig1,&nrow);
145 neig2 = read_xvg(Eig2File,&eig2,&nrow);
146 fprintf(stdout,"\nRead %d and %d eigenvalues\n",neig1,neig2);
147 n = min(n,min(neig1,neig2));
148 fprintf(stdout,"Will compare the covariance matrices using %d dimensions\n",
149 n);
150
151 sum1 = 0;
152 for(i=0; i<n; i++) {
153 if (eig1[1][i] < 0)
154 eig1[1][i] = 0;
155 sum1 += eig1[1][i];
156 eig1[1][i] = sqrt(eig1[1][i]);
157 }
158 trace1 = sum1;
159 for(i=n; i<neig1; i++)
160 trace1 += eig1[1][i];
161 sum2 = 0;
162 for(i=0; i<n; i++) {
163 if (eig2[1][i] < 0)
164 eig2[1][i] = 0;
165 sum2 += eig2[1][i];
166 eig2[1][i] = sqrt(eig2[1][i]);
167 }
168 trace2 = sum2;
169 for(i=n; i<neig2; i++)
170 trace2 += eig1[1][i];
171
172 fprintf(stdout,"Trace of the two matrices: %g and %g\n",sum1,sum2);
173 if (neig1!=n || neig2!=n)
174 fprintf(stdout,"this is %d%% and %d%% of the total trace\n",
175 (int)(100*sum1/trace1+0.5),(int)(100*sum2/trace2+0.5));
176 fprintf(stdout,"Square root of the traces: %g and %g\n",
177 sqrt(sum1),sqrt(sum2));
178
179 sab = 0;
180 for(i=0; i<n; i++) {
181 tmp = 0;
182 for(j=0; j<n; j++) {
183 ip = 0;
184 for(k=0; k<natoms; k++)
185 ip += iprod(eigvec1[i][k],eigvec2[j][k]);
186 tmp += eig2[1][j]*ip*ip;
187 }
188 sab += eig1[1][i]*tmp;
189 }
190
191 samsb2 = sum1+sum2-2*sab;
192 if (samsb2 < 0)
193 samsb2 = 0;
194
195 fprintf(stdout,"The overlap of the covariance matrices:\n");
196 fprintf(stdout," normalized: %.3f\n",1-sqrt(samsb2/(sum1+sum2)));
197 tmp = 1-sab/sqrt(sum1*sum2);
198 if (tmp < 0)
199 tmp = 0;
200 fprintf(stdout," shape: %.3f\n",1-sqrt(tmp));
201 }
202
203 static void inprod_matrix(char *matfile,int natoms,
204 int nvec1,int *eignr1,rvec **eigvec1,
205 int nvec2,int *eignr2,rvec **eigvec2,
206 bool bSelect,int noutvec,int *outvec)
207 {
208 FILE *out;
209 real **mat;
210 int i,x1,y1,x,y,nlevels;
211 int nx,ny;
212 real inp,*t_x,*t_y,max;
213 t_rgb rlo,rhi;
214
215 snew(t_y,nvec2);
216 if (bSelect) {
217 nx = noutvec;
218 ny = 0;
219 for(y1=0; y1<nx; y1++)
220 if (outvec[y1] < nvec2) {
221 t_y[ny] = eignr2[outvec[y1]]+1;
222 ny++;
223 }
224 } else {
225 nx = nvec1;
226 ny = nvec2;
227 for(y=0; y<ny; y++)
228 t_y[y] = eignr2[y]+1;
229 }
230
231 fprintf(stderr,"Calculating inner-product matrix of %dx%d eigenvectors\n",
232 nx,nvec2);
233
234 snew(mat,nx);
235 snew(t_x,nx);
236 max = 0;
237 for(x1=0; x1<nx; x1++) {
238 snew(mat[x1],ny);
239 if (bSelect)
240 x = outvec[x1];
241 else
242 x = x1;
243 t_x[x1] = eignr1[x]+1;
244 fprintf(stderr," %d",eignr1[x]+1);
245 for(y1=0; y1<ny; y1++) {
246 inp = 0;
247 if (bSelect) {
248 while (outvec[y1] >= nvec2)
249 y1++;
250 y= outvec[y1];
251 } else
252 y = y1;
253 for(i=0; i<natoms; i++)
254 inp += iprod(eigvec1[x][i],eigvec2[y][i]);
255 mat[x1][y1] = fabs(inp);
256 if (mat[x1][y1]>max)
257 max = mat[x1][y1];
258 }
259 }
260 fprintf(stderr,"\n");
261 rlo.r = 1; rlo.g = 1; rlo.b = 1;
262 rhi.r = 0; rhi.g = 0; rhi.b = 0;
263 nlevels = 41;
264 out = ffopen(matfile,"w");
265 write_xpm(out,"Eigenvector inner-products","in.prod.","run 1","run 2",
266 nx,ny,t_x,t_y,mat,0.0,max,rlo,rhi,&nlevels);
267 fclose(out);
268 }
269
270 static void overlap(char *outfile,int natoms,
271 rvec **eigvec1,
272 int nvec2,int *eignr2,rvec **eigvec2,
273 int noutvec,int *outvec)
274 {
275 FILE *out;
276 int i,v,vec,x;
277 real overlap,inp;
278
279 fprintf(stderr,"Calculating overlap between eigenvectors of set 2 with eigenvectors\n");
280 for(i=0; i<noutvec; i++)
281 fprintf(stderr,"%d ",outvec[i]+1);
282 fprintf(stderr,"\n");
283
284 out=xvgropen(outfile,"Subspace overlap",
285 "Eigenvectors of trajectory 2","Overlap");
286 fprintf(out,"@ subtitle \"using %d eigenvectors of trajectory 1\"\n",
287 noutvec);
288 overlap=0;
289 for(x=0; x<nvec2; x++) {
290 for(v=0; v<noutvec; v++) {
291 vec=outvec[v];
292 inp=0;
293 for(i=0; i<natoms; i++)
294 inp+=iprod(eigvec1[vec][i],eigvec2[x][i]);
295 overlap+=sqr(inp);
296 }
297 fprintf(out,"%5d %5.3f\n",eignr2[x]+1,overlap/noutvec);
298 }
299
300 fclose(out);
301 }
302
303 static void project(char *trajfile,t_topology *top,matrix topbox,rvec *xtop,
304 char *projfile,char *twodplotfile,char *threedplotfile,
305 char *filterfile,int skip,
306 char *extremefile,bool bExtrAll,real extreme,int nextr,
307 t_atoms *atoms,int natoms,atom_id *index,
308 bool bFit,rvec *xref,int nfit,atom_id *ifit,real *w_rls,
309 real *sqrtm,rvec *xav,
310 int *eignr,rvec **eigvec,
311 int noutvec,int *outvec, bool bSplit)
312 {
313 FILE *xvgrout=NULL;
314 int status,out=0,nat,i,j,d,v,vec,nfr,nframes=0,snew_size,frame;
315 int noutvec_extr,*imin,*imax;
316 atom_id *all_at;
317 matrix box;
318 rvec *xread,*x;
319 real t,inp,**inprod=NULL,min=0,max=0;
320 char str[STRLEN],str2[STRLEN],**ylabel,*c;
321
322 snew(x,natoms);
323
324 if (bExtrAll)
325 noutvec_extr=noutvec;
326 else
327 noutvec_extr=1;
328
329
330 if (trajfile) {
331 snew(inprod,noutvec+1);
332
333 if (filterfile) {
334 fprintf(stderr,"Writing a filtered trajectory to %s using eigenvectors\n",
335 filterfile);
336 for(i=0; i<noutvec; i++)
337 fprintf(stderr,"%d ",outvec[i]+1);
338 fprintf(stderr,"\n");
339 out=open_trx(filterfile,"w");
340 }
341 snew_size=0;
342 nfr=0;
343 nframes=0;
344 nat=read_first_x(&status,trajfile,&t,&xread,box);
345 if (nat>atoms->nr)
346 fatal_error(0,"the number of atoms in your trajectory (%d) is larger than the number of atoms in your structure file (%d)",nat,atoms->nr);
347 snew(all_at,nat);
348 for(i=0; i<nat; i++)
349 all_at[i]=i;
350 do {
351 if (nfr % skip == 0) {
352 if (top)
353 rm_pbc(&(top->idef),nat,box,xread,xread);
354 if (nframes>=snew_size) {
355 snew_size+=100;
356 for(i=0; i<noutvec+1; i++)
357 srenew(inprod[i],snew_size);
358 }
359 inprod[noutvec][nframes]=t;
360 /* calculate x: a fitted struture of the selected atoms */
361 if (bFit && (xref==NULL)) {
362 reset_x(nfit,ifit,nat,NULL,xread,w_rls);
363 do_fit(nat,w_rls,xtop,xread);
364 }
365 for (i=0; i<natoms; i++)
366 copy_rvec(xread[index[i]],x[i]);
367 if (bFit && xref) {
368 reset_x(natoms,all_at,natoms,NULL,x,w_rls);
369 do_fit(natoms,w_rls,xref,x);
370 }
371
372 for(v=0; v<noutvec; v++) {
373 vec=outvec[v];
374 /* calculate (mass-weighted) projection */
375 inp=0;
376 for (i=0; i<natoms; i++) {
377 inp+=(eigvec[vec][i][0]*(x[i][0]-xav[i][0])+
378 eigvec[vec][i][1]*(x[i][1]-xav[i][1])+
379 eigvec[vec][i][2]*(x[i][2]-xav[i][2]))*sqrtm[i];
380 }
381 inprod[v][nframes]=inp;
382 }
383 if (filterfile) {
384 for(i=0; i<natoms; i++)
385 for(d=0; d<DIM; d++) {
386 /* misuse xread for output */
387 xread[index[i]][d] = xav[i][d];
388 for(v=0; v<noutvec; v++)
389 xread[index[i]][d] +=
390 inprod[v][nframes]*eigvec[outvec[v]][i][d]/sqrtm[i];
391 }
392 write_trx(out,natoms,index,atoms,0,t,box,xread,NULL);
393 }
394 nframes++;
395 }
396 nfr++;
397 } while (read_next_x(status,&t,nat,xread,box));
398 close_trj(status);
399 sfree(x);
400 if (filterfile)
401 close_trx(out);
402 }
403 else
404 snew(xread,atoms->nr);
405
406
407 if (projfile) {
408 snew(ylabel,noutvec);
409 for(v=0; v<noutvec; v++) {
410 sprintf(str,"vec %d",eignr[outvec[v]]+1);
411 ylabel[v]=strdup(str);
412 }
413 sprintf(str,"projection on eigenvectors (%s)",proj_unit);
414 write_xvgr_graphs(projfile, noutvec, str, xvgr_tlabel(),
415 ylabel, nframes, inprod[noutvec], inprod,
416 time_factor(), FALSE, bSplit);
417 }
418
419 if (twodplotfile) {
420 sprintf(str,"projection on eigenvector %d (%s)",
421 eignr[outvec[0]]+1,proj_unit);
422 sprintf(str2,"projection on eigenvector %d (%s)",
423 eignr[outvec[noutvec-1]]+1,proj_unit);
424 xvgrout=xvgropen(twodplotfile,"2D projection of trajectory",str,str2);
425 for(i=0; i<nframes; i++) {
426 if ( bSplit && i>0 && abs(inprod[noutvec][i])<1e-5 )
427 fprintf(xvgrout,"&\n");
428 fprintf(xvgrout,"%10.5f %10.5f\n",inprod[0][i],inprod[noutvec-1][i]);
429 }
430 fclose(xvgrout);
431 }
432
433 if (threedplotfile) {
434 t_atoms atoms;
435 rvec *x;
436 real *b=NULL;
437 matrix box;
438 char *resnm,*atnm, pdbform[STRLEN];
439 bool bPDB, b4D;
440 FILE *out;
441
442 if (noutvec < 3)
443 fatal_error(0,"You have selected less than 3 eigenvectors");
444
445 /* initialize */
446 bPDB = fn2ftp(threedplotfile)==efPDB;
447 clear_mat(box);
448 box[XX][XX] = box[YY][YY] = box[ZZ][ZZ] = 1;
449
450 b4D = bPDB && (noutvec >= 4);
451 if (b4D) {
452 fprintf(stderr, "You have selected four or more eigenvectors:\n"
453 "fourth eigenvector will be plotted "
454 "in bfactor field of pdb file\n");
455 sprintf(str,"4D proj. of traj. on eigenv. %d, %d, %d and %d",
456 eignr[outvec[0]]+1,eignr[outvec[1]]+1,
457 eignr[outvec[2]]+1,eignr[outvec[3]]+1);
458 } else {
459 sprintf(str,"3D proj. of traj. on eigenv. %d, %d and %d",
460 eignr[outvec[0]]+1,eignr[outvec[1]]+1,eignr[outvec[2]]+1);
461 }
462 init_t_atoms(&atoms,nframes,FALSE);
463 snew(x,nframes);
464 snew(b,nframes);
465 atnm=strdup("C");
466 resnm=strdup("PRJ");
467 for(i=0; i<nframes; i++) {
468 atoms.resname[i]=&resnm;
469 atoms.atomname[i]=&atnm;
470 atoms.atom[i].resnr=i;
471 x[i][XX]=inprod[0][i];
472 x[i][YY]=inprod[1][i];
473 x[i][ZZ]=inprod[2][i];
474 if (b4D)
475 b[i] =inprod[3][i];
476 }
477 if ( ( b4D || bSplit ) && bPDB ) {
478 strcpy(pdbform,pdbformat);
479 strcat(pdbform,"%8.4f%8.4f\n");
480
481 out=ffopen(threedplotfile,"w");
482 fprintf(out,"HEADER %s\n",str);
483 if ( b4D )
484 fprintf(out,"REMARK %s\n","fourth dimension plotted as B-factor");
485 j=0;
486 for(i=0; i<atoms.nr; i++) {
487 if ( j>0 && bSplit && abs(inprod[noutvec][i])<1e-5 ) {
488 fprintf(out,"TER\n");
489 j=0;
490 }
491 fprintf(out,pdbform,"ATOM",i+1,"C","PRJ",' ',j+1,
492 PR_VEC(10*x[i]), 1.0, 10*b[i]);
493 if (j>0)
494 fprintf(out,"CONECT%5d%5d\n", i, i+1);
495 j++;
496 }
497 fprintf(out,"TER\n");
498 fclose(out);
499 } else
500 write_sto_conf(threedplotfile,str,&atoms,x,NULL,box);
501 free_t_atoms(&atoms);
502 }
503
504 if (extremefile) {
505 if (extreme==0) {
506 fprintf(stderr,"%11s %17s %17s\n","eigenvector","Minimum","Maximum");
507 fprintf(stderr,
508 "%11s %10s %10s %10s %10s\n","","value","time","value","time");
509 snew(imin,noutvec_extr);
510 snew(imax,noutvec_extr);
511 for(v=0; v<noutvec_extr; v++) {
512 for(i=0; i<nframes; i++) {
513 if (inprod[v][i]<inprod[v][imin[v]])
514 imin[v]=i;
515 if (inprod[v][i]>inprod[v][imax[v]])
516 imax[v]=i;
517 }
518 min=inprod[v][imin[v]];
519 max=inprod[v][imax[v]];
520 fprintf(stderr,"%7d %10.6f %10.1f %10.6f %10.1f\n",
521 eignr[outvec[v]]+1,
522 min,inprod[noutvec][imin[v]],max,inprod[noutvec][imax[v]]);
523 }
524 }
525 else {
526 min=-extreme;
527 max=+extreme;
528 }
529 /* build format string for filename: */
530 strcpy(str,extremefile);/* copy filename */
531 c=strrchr(str,'.'); /* find where extention begins */
532 strcpy(str2,c); /* get extention */
533 sprintf(c,"%%d%s",str2); /* append '%s' and extention to filename */
534 for(v=0; v<noutvec_extr; v++) {
535 /* make filename using format string */
536 if (noutvec_extr==1)
537 strcpy(str2,extremefile);
538 else
539 sprintf(str2,str,eignr[outvec[v]]+1);
540 fprintf(stderr,"Writing %d frames along eigenvector %d to %s\n",
541 nextr,outvec[v]+1,str2);
542 out=open_trx(str2,"w");
543 for(frame=0; frame<nextr; frame++) {
544 if ((extreme==0) && (nextr<=3))
545 for(i=0; i<natoms; i++)
546 atoms->atom[index[i]].chain='A'+frame;
547 for(i=0; i<natoms; i++)
548 for(d=0; d<DIM; d++)
549 xread[index[i]][d] =
550 (xav[i][d] + (min*(nextr-frame-1)+max*frame)/(nextr-1)
551 *eigvec[outvec[v]][i][d]/sqrtm[i]);
552 write_trx(out,natoms,index,atoms,0,frame,topbox,xread,NULL);
553 }
554 close_trx(out);
555 }
556 }
557 fprintf(stderr,"\n");
558 }
559
560 static void components(char *outfile,int natoms,real *sqrtm,
561 int *eignr,rvec **eigvec,
562 int noutvec,int *outvec)
563 {
564 int g,v,i;
565 real *x,**y;
566 char str[STRLEN],**ylabel;
567
568 fprintf(stderr,"Writing atom displacements to %s\n",outfile);
569
570 snew(ylabel,noutvec);
571 snew(y,noutvec);
572 snew(x,natoms);
573 for(i=0; i<natoms; i++)
574 x[i]=i+1;
575 for(g=0; g<noutvec; g++) {
576 v=outvec[g];
577 sprintf(str,"vec %d",eignr[v]+1);
578 ylabel[g]=strdup(str);
579 snew(y[g],natoms);
580 for(i=0; i<natoms; i++)
581 y[g][i] = norm(eigvec[v][i])/sqrtm[i];
582 }
583 write_xvgr_graphs(outfile,noutvec,"Atom displacements (nm)","Atom number",
584 ylabel,natoms,x,y,1,TRUE,FALSE);
585 fprintf(stderr,"\n");
586 }
587
588 int main(int argc,char *argv[])
589 {
590 static char *desc[] = {
591 "[TT]g_anaeig[tt] analyzes eigenvectors. The eigenvectors can be of a",
592 "covariance matrix ([TT]g_covar[tt]) or of a Normal Modes anaysis",
593 "([TT]g_nmeig[tt]).[PAR]",
594
595 "When a trajectory is projected on eigenvectors, all structures are",
596 "fitted to the structure in the eigenvector file, if present, otherwise",
597 "to the structure in the structure file. When no run input file is",
598 "supplied, periodicity will not be taken into account. Most analyses",
599 "are performed on eigenvectors [TT]-first[tt] to [TT]-last[tt], but when",
600 "[TT]-first[tt] is set to -1 you will be prompted for a selection.[PAR]",
601
602 "[TT]-disp[tt]: plot all atom displacements of eigenvectors",
603 "[TT]-first[tt] to [TT]-last[tt].[PAR]",
604
605 "[TT]-proj[tt]: calculate projections of a trajectory on eigenvectors",
606 "[TT]-first[tt] to [TT]-last[tt].",
607 "The projections of a trajectory on the eigenvectors of its",
608 "covariance matrix are called principal components (pc's).",
609 "It is often useful to check the cosine content the pc's,",
610 "since the pc's of random diffusion are cosines with the number",
611 "of periods equal to half the pc index.",
612 "The cosine content of the pc's can be calculated with the program",
613 "[TT]g_analyze[tt].[PAR]",
614
615 "[TT]-2d[tt]: calculate a 2d projection of a trajectory on eigenvectors",
616 "[TT]-first[tt] and [TT]-last[tt].[PAR]",
617
618 "[TT]-3d[tt]: calculate a 3d projection of a trajectory on the first",
619 "three selected eigenvectors.[PAR]",
620
621 "[TT]-filt[tt]: filter the trajectory to show only the motion along",
622 "eigenvectors [TT]-first[tt] to [TT]-last[tt].[PAR]",
623
624 "[TT]-extr[tt]: calculate the two extreme projections along a trajectory",
625 "on the average structure and interpolate [TT]-nframes[tt] frames",
626 "between them, or set your own extremes with [TT]-max[tt]. The",
627 "eigenvector [TT]-first[tt] will be written unless [TT]-first[tt] and",
628 "[TT]-last[tt] have been set explicitly, in which case all eigenvectors",
629 "will be written to separate files. Chain identifiers will be added",
630 "when writing a [TT].pdb[tt] file with two or three structures (you",
631 "can use [TT]rasmol -nmrpdb[tt] to view such a pdb file).[PAR]",
632
633 " Overlap calculations between covariance analysis:[BR]",
634 " NOTE: the analysis should use the same fitting structure[PAR]",
635
636 "[TT]-over[tt]: calculate the subspace overlap of the eigenvectors in",
637 "file [TT]-v2[tt] with eigenvectors [TT]-first[tt] to [TT]-last[tt]",
638 "in file [TT]-v[tt].[PAR]",
639
640 "[TT]-inpr[tt]: calculate a matrix of inner-products between",
641 "eigenvectors in files [TT]-v[tt] and [TT]-v2[tt]. All eigenvectors",
642 "of both files will be used unless [TT]-first[tt] and [TT]-last[tt]",
643 "have been set explicitly.[PAR]",
644
645 "When [TT]-v[tt], [TT]-eig1[tt], [TT]-v2[tt] and [TT]-eig2[tt] are given,",
646 "a single number for the overlap between the covariance matrices is",
647 "generated. The formulas are:[BR]",
648 " difference = sqrt(tr((sqrt(M1) - sqrt(M2))^2))[BR]",
649 "normalized overlap = 1 - difference/sqrt(tr(M1) + tr(M2))[BR]",
650 " shape overlap = 1 - sqrt(tr((sqrt(M1/tr(M1)) - sqrt(M2/tr(M2)))^2))[BR]",
651 "where M1 and M2 are the two covariance matrices and tr is the trace",
652 "of a matrix. The numbers are proportional to the overlap of the square",
653 "root of the fluctuations. The normalized overlap is the most useful",
654 "number, it is 1 for identical matrices and 0 when the sampled",
655 "subspaces are orthogonal."
656 };
657 static int first=1,last=8,skip=1,nextr=2;
658 static real max=0.0;
659 static bool bSplit=FALSE;
660 t_pargs pa[] = {
661 { "-first", FALSE, etINT, {&first},
662 "First eigenvector for analysis (-1 is select)" },
663 { "-last", FALSE, etINT, {&last},
664 "Last eigenvector for analysis (-1 is till the last)" },
665 { "-skip", FALSE, etINT, {&skip},
666 "Only analyse every nr-th frame" },
667 { "-max", FALSE, etREAL, {&max},
668 "Maximum for projection of the eigenvector on the average structure, "
669 "max=0 gives the extremes" },
670 { "-nframes", FALSE, etINT, {&nextr},
671 "Number of frames for the extremes output" },
672 { "-split", FALSE, etBOOL, {&bSplit},
673 "Split eigenvector projections where time is zero" }
674 };
675 #define NPA asize(pa)
676
677 FILE *out;
678 int status,trjout;
679 t_topology top;
680 t_atoms *atoms=NULL;
681 rvec *xtop,*xref1,*xref2;
682 bool bDMR1,bDMA1,bDMR2,bDMA2;
683 int nvec1,nvec2,*eignr1=NULL,*eignr2=NULL;
684 rvec *x,*xread,*xav1,*xav2,**eigvec1=NULL,**eigvec2=NULL;
685 matrix topbox;
686 real xid,totmass,*sqrtm,*w_rls,t,lambda;
687 int natoms,step;
688 char *grpname,*indexfile,title[STRLEN];
689 int i,j,d;
690 int nout,*iout,noutvec,*outvec,nfit;
691 atom_id *index,*ifit;
692 char *Vec2File,*topfile;
693 char *Eig1File,*Eig2File;
694 char *CompFile,*ProjOnVecFile;
695 char *TwoDPlotFile,*ThreeDPlotFile;
696 char *FilterFile,*ExtremeFile;
697 char *OverlapFile,*InpMatFile;
698 bool bFit1,bFit2,bM,bIndex,bTPS,bTop,bVec2,bProj;
699 bool bFirstToLast,bFirstLastSet,bTraj,bCompare,bPDB3D;
700
701 t_filenm fnm[] = {
702 { efTRN, "-v", "eigenvec", ffREAD },
703 { efTRN, "-v2", "eigenvec2", ffOPTRD },
704 { efTRX, "-f", NULL, ffOPTRD },
705 { efTPS, NULL, NULL, ffOPTRD },
706 { efNDX, NULL, NULL, ffOPTRD },
707 { efXVG, "-eig1", "eigenval1", ffOPTRD },
708 { efXVG, "-eig2", "eigenval2", ffOPTRD },
709 { efXVG, "-disp", "eigdisp", ffOPTWR },
710 { efXVG, "-proj", "proj", ffOPTWR },
711 { efXVG, "-2d", "2dproj", ffOPTWR },
712 { efSTO, "-3d", "3dproj.pdb", ffOPTWR },
713 { efTRX, "-filt", "filtered", ffOPTWR },
714 { efTRX, "-extr", "extreme.pdb", ffOPTWR },
715 { efXVG, "-over", "overlap", ffOPTWR },
716 { efXPM, "-inpr", "inprod", ffOPTWR }
717 };
718 #define NFILE asize(fnm)
719
720 CopyRight(stderr,argv[0]);
721 parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW | PCA_BE_NICE ,
722 NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL);
723
724 indexfile=ftp2fn_null(efNDX,NFILE,fnm);
725
726 Vec2File = opt2fn_null("-v2",NFILE,fnm);
727 topfile = ftp2fn(efTPS,NFILE,fnm);
728 Eig1File = opt2fn_null("-eig1",NFILE,fnm);
729 Eig2File = opt2fn_null("-eig2",NFILE,fnm);
730 CompFile = opt2fn_null("-disp",NFILE,fnm);
731 ProjOnVecFile = opt2fn_null("-proj",NFILE,fnm);
732 TwoDPlotFile = opt2fn_null("-2d",NFILE,fnm);
733 ThreeDPlotFile = opt2fn_null("-3d",NFILE,fnm);
734 FilterFile = opt2fn_null("-filt",NFILE,fnm);
735 ExtremeFile = opt2fn_null("-extr",NFILE,fnm);
736 OverlapFile = opt2fn_null("-over",NFILE,fnm);
737 InpMatFile = ftp2fn_null(efXPM,NFILE,fnm);
738 bTop = fn2bTPX(topfile);
739 bProj = ProjOnVecFile || TwoDPlotFile || ThreeDPlotFile
740 || FilterFile || ExtremeFile;
741 bFirstLastSet =
742 opt2parg_bSet("-first",NPA,pa) && opt2parg_bSet("-last",NPA,pa);
743 bFirstToLast = CompFile || ProjOnVecFile || FilterFile || OverlapFile ||
744 ((ExtremeFile || InpMatFile) && bFirstLastSet);
745 bVec2 = Vec2File || OverlapFile || InpMatFile;
746 bM = CompFile || bProj;
747 bTraj = ProjOnVecFile || FilterFile || (ExtremeFile && (max==0))
748 || TwoDPlotFile || ThreeDPlotFile;
749 bIndex = bM || bProj;
750 bTPS = ftp2bSet(efTPS,NFILE,fnm) || bM || bTraj ||
751 FilterFile || (bIndex && indexfile);
752 bCompare = Vec2File && Eig1File && Eig2File;
753 bPDB3D = fn2ftp(ThreeDPlotFile)==efPDB;
754
755 read_eigenvectors(opt2fn("-v",NFILE,fnm),&natoms,&bFit1,
756 &xref1,&bDMR1,&xav1,&bDMA1,&nvec1,&eignr1,&eigvec1);
757 if (bVec2) {
758 read_eigenvectors(Vec2File,&i,&bFit2,
759 &xref2,&bDMR2,&xav2,&bDMA2,&nvec2,&eignr2,&eigvec2);
760 if (i!=natoms)
761 fatal_error(0,"Dimensions in the eigenvector files don't match");
762 }
763
764 if ((!bFit1 || xref1) && !bDMR1 && !bDMA1)
765 bM=FALSE;
766 if ((xref1==NULL) && (bM || bTraj))
767 bTPS=TRUE;
768
769 xtop=NULL;
770 nfit=0;
771 ifit=NULL;
772 w_rls=NULL;
773 if (!bTPS)
774 bTop=FALSE;
775 else {
776 bTop=read_tps_conf(ftp2fn(efTPS,NFILE,fnm),
777 title,&top,&xtop,NULL,topbox,bM);
778 atoms=&top.atoms;
779 rm_pbc(&(top.idef),atoms->nr,topbox,xtop,xtop);
780 /* Fitting is only needed when we need to read a trajectory */
781 if (bTraj) {
782 if ((xref1==NULL) || (bM && bDMR1)) {
783 if (bFit1) {
784 printf("\nNote: the structure in %s should be the same\n"
785 " as the one used for the fit in g_covar\n",topfile);
786 printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
787 get_index(atoms,indexfile,1,&nfit,&ifit,&grpname);
788 snew(w_rls,atoms->nr);
789 for(i=0; (i<nfit); i++)
790 if (bM && bDMR1)
791 w_rls[ifit[i]]=atoms->atom[ifit[i]].m;
792 else
793 w_rls[ifit[i]]=1.0;
794 }
795 else {
796 /* make the fit index in xref instead of xtop */
797 nfit=natoms;
798 snew(ifit,natoms);
799 snew(w_rls,nfit);
800 for(i=0; (i<nfit); i++) {
801 ifit[i]=i;
802 w_rls[i]=atoms->atom[ifit[i]].m;
803 }
804 }
805 }
806 else {
807 /* make the fit non mass weighted on xref */
808 nfit=natoms;
809 snew(ifit,nfit);
810 snew(w_rls,nfit);
811 for(i=0; i<nfit; i++) {
812 ifit[i]=i;
813 w_rls[i]=1.0;
814 }
815 }
816 }
817 }
818
819 if (bIndex) {
820 printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n",natoms);
821 get_index(atoms,indexfile,1,&i,&index,&grpname);
822 if (i!=natoms)
823 fatal_error(0,"you selected a group with %d elements instead of %d",
824 i,natoms);
825 printf("\n");
826 }
827
828 snew(sqrtm,natoms);
829 if (bM && bDMA1) {
830 proj_unit="u\\S1/2\\Nnm";
831 for(i=0; (i<natoms); i++)
832 sqrtm[i]=sqrt(atoms->atom[index[i]].m);
833 } else {
834 proj_unit="nm";
835 for(i=0; (i<natoms); i++)
836 sqrtm[i]=1.0;
837 }
838
839 if (bVec2) {
840 t=0;
841 totmass=0;
842 for(i=0; (i<natoms); i++)
843 for(d=0;(d<DIM);d++) {
844 t+=sqr((xav1[i][d]-xav2[i][d])*sqrtm[i]);
845 totmass+=sqr(sqrtm[i]);
846 }
847 fprintf(stdout,"RMSD (without fit) between the two average structures:"
848 " %.3f (nm)\n\n",sqrt(t/totmass));
849 }
850
851 if (last==-1)
852 last=natoms*DIM;
853 if (first>-1) {
854 if (bFirstToLast) {
855 /* make an index from first to last */
856 nout=last-first+1;
857 snew(iout,nout);
858 for(i=0; i<nout; i++)
859 iout[i]=first-1+i;
860 } else if (ThreeDPlotFile) {
861 /* make an index of first+(0,1,2) and last */
862 nout = bPDB3D ? 4 : 3;
863 nout = min(last-first+1, nout);
864 snew(iout,nout);
865 iout[0]=first-1;
866 iout[1]=first;
867 if (nout>3)
868 iout[2]=first+1;
869 iout[nout-1]=last-1;
870 } else {
871 /* make an index of first and last */
872 nout=2;
873 snew(iout,nout);
874 iout[0]=first-1;
875 iout[1]=last-1;
876 }
877 }
878 else {
879 printf("Select eigenvectors for output, end your selection with 0\n");
880 nout=-1;
881 iout=NULL;
882 do {
883 nout++;
884 srenew(iout,nout+1);
885 scanf("%d",&iout[nout]);
886 iout[nout]--;
887 } while (iout[nout]>=0);
888 printf("\n");
889 }
890 /* make an index of the eigenvectors which are present */
891 snew(outvec,nout);
892 noutvec=0;
893 for(i=0; i<nout; i++) {
894 j=0;
895 while ((j<nvec1) && (eignr1[j]!=iout[i]))
896 j++;
897 if ((j<nvec1) && (eignr1[j]==iout[i])) {
898 outvec[noutvec]=j;
899 noutvec++;
900 }
901 }
902 fprintf(stderr,"%d eigenvectors selected for output",noutvec);
903 if (noutvec <= 100) {
904 fprintf(stderr,":");
905 for(j=0; j<noutvec; j++)
906 fprintf(stderr," %d",eignr1[outvec[j]]+1);
907 }
908 fprintf(stderr,"\n");
909
910 if (CompFile)
911 components(CompFile,natoms,sqrtm,eignr1,eigvec1,noutvec,outvec);
912
913 if (bProj)
914 project(bTraj ? opt2fn("-f",NFILE,fnm) : NULL,
915 bTop ? &top : NULL,topbox,xtop,
916 ProjOnVecFile,TwoDPlotFile,ThreeDPlotFile,FilterFile,skip,
917 ExtremeFile,bFirstLastSet,max,nextr,atoms,natoms,index,
918 bFit1,xref1,nfit,ifit,w_rls,
919 sqrtm,xav1,eignr1,eigvec1,noutvec,outvec,bSplit);
920
921 if (OverlapFile)
922 overlap(OverlapFile,natoms,
923 eigvec1,nvec2,eignr2,eigvec2,noutvec,outvec);
924
925 if (InpMatFile)
926 inprod_matrix(InpMatFile,natoms,
927 nvec1,eignr1,eigvec1,nvec2,eignr2,eigvec2,
928 bFirstLastSet,noutvec,outvec);
929
930 if (bCompare)
931 compare(natoms,nvec1,eigvec1,nvec2,eigvec2,Eig1File,Eig2File);
932
933 if (!CompFile && !bProj && !OverlapFile && !InpMatFile && !bCompare)
934 fprintf(stderr,"\nIf you want some output,"
935 " set one (or two or ...) of the output file options\n");
936
937 view_all(NFILE, fnm);
938
939 return 0;
940 }
941
The ultimate molecular dynamics simulation package