3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
45 #include "gmx_fatal.h"
65 static void calc_entropy_qh(FILE *fp,int n,real eigval[],real temp,int nskip)
71 hbar = PLANCK1/(2*M_PI);
72 for(i=0; (i<n-nskip); i++) {
74 lambda = eigval[i]*AMU;
75 w = sqrt(BOLTZMANN*temp/lambda)/NANO;
76 hwkT = (hbar*w)/(BOLTZMANN*temp);
77 dS = (hwkT/(exp(hwkT) - 1) - log(1-exp(-hwkT)));
80 fprintf(debug,"i = %5d w = %10g lam = %10g hwkT = %10g dS = %10g\n",
84 fprintf(stderr,"eigval[%d] = %g\n",i,eigval[i]);
88 fprintf(fp,"The Entropy due to the Quasi Harmonic approximation is %g J/mol K\n",
92 static void calc_entropy_schlitter(FILE *fp,int n,int nskip,
93 real *eigval,real temp)
99 double hbar,kt,kteh,S;
101 hbar = PLANCK1/(2*M_PI);
103 kteh = kt*exp(2.0)/(hbar*hbar)*AMU*(NANO*NANO);
105 fprintf(debug,"n = %d, nskip = %d kteh = %g\n",n,nskip,kteh);
108 for(i=0; (i<n-nskip); i++) {
109 dd = 1+kteh*eigval[i];
114 fprintf(fp,"The Entropy due to the Schlitter formula is %g J/mol K\n",S);
117 const char *proj_unit;
119 static real tick_spacing(real range,int minticks)
127 sp = 0.2*exp(log(10)*ceil(log(range)/log(10)));
128 while (range/sp < minticks-1)
134 static void write_xvgr_graphs(const char *file, int ngraphs, int nsetspergraph,
135 const char *title, const char *subtitle,
136 const char *xlabel, const char **ylabel,
137 int n, real *x, real **y, real ***sy,
138 real scale_x, gmx_bool bZero, gmx_bool bSplit,
139 const output_env_t oenv)
143 real min,max,xsp,ysp;
145 out=ffopen(file,"w");
146 if (output_env_get_xvg_format(oenv) == exvgXMGRACE) {
147 fprintf(out,"@ autoscale onread none\n");
149 for(g=0; g<ngraphs; g++) {
154 if (y[g][i]<min) min=y[g][i];
155 if (y[g][i]>max) max=y[g][i];
160 for(s=0; s<nsetspergraph; s++)
162 if (sy[g][s][i]<min) min=sy[g][s][i];
163 if (sy[g][s][i]>max) max=sy[g][s][i];
169 min=min-0.1*(max-min);
170 max=max+0.1*(max-min);
171 xsp=tick_spacing((x[n-1]-x[0])*scale_x,4);
172 ysp=tick_spacing(max-min,3);
173 if (output_env_get_print_xvgr_codes(oenv)) {
174 fprintf(out,"@ with g%d\n@ g%d on\n",g,g);
176 fprintf(out,"@ title \"%s\"\n",title);
178 fprintf(out,"@ subtitle \"%s\"\n",subtitle);
181 fprintf(out,"@ xaxis label \"%s\"\n",xlabel);
183 fprintf(out,"@ xaxis ticklabel off\n");
185 fprintf(out,"@ world xmin %g\n",x[0]*scale_x);
186 fprintf(out,"@ world xmax %g\n",x[n-1]*scale_x);
187 fprintf(out,"@ world ymin %g\n",min);
188 fprintf(out,"@ world ymax %g\n",max);
190 fprintf(out,"@ view xmin 0.15\n");
191 fprintf(out,"@ view xmax 0.85\n");
192 fprintf(out,"@ view ymin %g\n",0.15+(ngraphs-1-g)*0.7/ngraphs);
193 fprintf(out,"@ view ymax %g\n",0.15+(ngraphs-g)*0.7/ngraphs);
194 fprintf(out,"@ yaxis label \"%s\"\n",ylabel[g]);
195 fprintf(out,"@ xaxis tick major %g\n",xsp);
196 fprintf(out,"@ xaxis tick minor %g\n",xsp/2);
197 fprintf(out,"@ xaxis ticklabel start type spec\n");
198 fprintf(out,"@ xaxis ticklabel start %g\n",ceil(min/xsp)*xsp);
199 fprintf(out,"@ yaxis tick major %g\n",ysp);
200 fprintf(out,"@ yaxis tick minor %g\n",ysp/2);
201 fprintf(out,"@ yaxis ticklabel start type spec\n");
202 fprintf(out,"@ yaxis ticklabel start %g\n",ceil(min/ysp)*ysp);
203 if ((min<0) && (max>0)) {
204 fprintf(out,"@ zeroxaxis bar on\n");
205 fprintf(out,"@ zeroxaxis bar linestyle 3\n");
208 for(s=0; s<nsetspergraph; s++) {
210 if ( bSplit && i>0 && abs(x[i])<1e-5 )
212 if (output_env_get_print_xvgr_codes(oenv))
215 fprintf(out,"%10.4f %10.5f\n",
216 x[i]*scale_x,y ? y[g][i] : sy[g][s][i]);
218 if (output_env_get_print_xvgr_codes(oenv))
226 compare(int natoms,int n1,rvec **eigvec1,int n2,rvec **eigvec2,
227 real *eigval1, int neig1, real *eigval2, int neig2)
231 double sum1,sum2,trace1,trace2,sab,samsb2,tmp,ip;
235 n = min(n,min(neig1,neig2));
236 fprintf(stdout,"Will compare the covariance matrices using %d dimensions\n",n);
244 eigval1[i] = sqrt(eigval1[i]);
247 for(i=n; i<neig1; i++)
248 trace1 += eigval1[i];
255 eigval2[i] = sqrt(eigval2[i]);
258 for(i=n; i<neig2; i++)
259 trace2 += eigval2[i];
261 fprintf(stdout,"Trace of the two matrices: %g and %g\n",sum1,sum2);
262 if (neig1!=n || neig2!=n)
263 fprintf(stdout,"this is %d%% and %d%% of the total trace\n",
264 (int)(100*sum1/trace1+0.5),(int)(100*sum2/trace2+0.5));
265 fprintf(stdout,"Square root of the traces: %g and %g\n",
266 sqrt(sum1),sqrt(sum2));
275 for(k=0; k<natoms; k++)
276 ip += iprod(eigvec1[i][k],eigvec2[j][k]);
277 tmp += eigval2[j]*ip*ip;
279 sab += eigval1[i]*tmp;
282 samsb2 = sum1+sum2-2*sab;
286 fprintf(stdout,"The overlap of the covariance matrices:\n");
287 fprintf(stdout," normalized: %.3f\n",1-sqrt(samsb2/(sum1+sum2)));
288 tmp = 1-sab/sqrt(sum1*sum2);
291 fprintf(stdout," shape: %.3f\n",1-sqrt(tmp));
295 static void inprod_matrix(const char *matfile,int natoms,
296 int nvec1,int *eignr1,rvec **eigvec1,
297 int nvec2,int *eignr2,rvec **eigvec2,
298 gmx_bool bSelect,int noutvec,int *outvec)
302 int i,x1,y1,x,y,nlevels;
304 real inp,*t_x,*t_y,max;
311 for(y1=0; y1<nx; y1++)
312 if (outvec[y1] < nvec2) {
313 t_y[ny] = eignr2[outvec[y1]]+1;
320 t_y[y] = eignr2[y]+1;
323 fprintf(stderr,"Calculating inner-product matrix of %dx%d eigenvectors\n",
329 for(x1=0; x1<nx; x1++) {
335 t_x[x1] = eignr1[x]+1;
336 fprintf(stderr," %d",eignr1[x]+1);
337 for(y1=0; y1<ny; y1++) {
340 while (outvec[y1] >= nvec2)
345 for(i=0; i<natoms; i++)
346 inp += iprod(eigvec1[x][i],eigvec2[y][i]);
347 mat[x1][y1] = fabs(inp);
352 fprintf(stderr,"\n");
353 rlo.r = 1; rlo.g = 1; rlo.b = 1;
354 rhi.r = 0; rhi.g = 0; rhi.b = 0;
356 out = ffopen(matfile,"w");
357 write_xpm(out,0,"Eigenvector inner-products","in.prod.","run 1","run 2",
358 nx,ny,t_x,t_y,mat,0.0,max,rlo,rhi,&nlevels);
362 static void overlap(const char *outfile,int natoms,
364 int nvec2,int *eignr2,rvec **eigvec2,
365 int noutvec,int *outvec,
366 const output_env_t oenv)
372 fprintf(stderr,"Calculating overlap between eigenvectors of set 2 with eigenvectors\n");
373 for(i=0; i<noutvec; i++)
374 fprintf(stderr,"%d ",outvec[i]+1);
375 fprintf(stderr,"\n");
377 out=xvgropen(outfile,"Subspace overlap",
378 "Eigenvectors of trajectory 2","Overlap",oenv);
379 fprintf(out,"@ subtitle \"using %d eigenvectors of trajectory 1\"\n",
382 for(x=0; x<nvec2; x++) {
383 for(v=0; v<noutvec; v++) {
386 for(i=0; i<natoms; i++)
387 inp+=iprod(eigvec1[vec][i],eigvec2[x][i]);
390 fprintf(out,"%5d %5.3f\n",eignr2[x]+1,overlap/noutvec);
396 static void project(const char *trajfile,t_topology *top,int ePBC,matrix topbox,
397 const char *projfile,const char *twodplotfile,
398 const char *threedplotfile, const char *filterfile,int skip,
399 const char *extremefile,gmx_bool bExtrAll,real extreme,
400 int nextr, t_atoms *atoms,int natoms,atom_id *index,
401 gmx_bool bFit,rvec *xref,int nfit,atom_id *ifit,real *w_rls,
402 real *sqrtm,rvec *xav,
403 int *eignr,rvec **eigvec,
404 int noutvec,int *outvec, gmx_bool bSplit,
405 const output_env_t oenv)
408 int nat,i,j,d,v,vec,nfr,nframes=0,snew_size,frame;
409 t_trxstatus *out=NULL;
411 int noutvec_extr,*imin,*imax;
415 real t,inp,**inprod=NULL,min=0,max=0;
416 char str[STRLEN],str2[STRLEN],**ylabel,*c;
418 gmx_rmpbc_t gpbc=NULL;
423 noutvec_extr=noutvec;
429 snew(inprod,noutvec+1);
432 fprintf(stderr,"Writing a filtered trajectory to %s using eigenvectors\n",
434 for(i=0; i<noutvec; i++)
435 fprintf(stderr,"%d ",outvec[i]+1);
436 fprintf(stderr,"\n");
437 out=open_trx(filterfile,"w");
442 nat=read_first_x(oenv,&status,trajfile,&t,&xread,box);
444 gmx_fatal(FARGS,"the number of atoms in your trajectory (%d) is larger than the number of atoms in your structure file (%d)",nat,atoms->nr);
448 gpbc = gmx_rmpbc_init(&top->idef,ePBC,nat,box);
450 gmx_rmpbc_done(gpbc);
455 if (nfr % skip == 0) {
457 gmx_rmpbc(gpbc,nat,box,xread);
458 if (nframes>=snew_size) {
460 for(i=0; i<noutvec+1; i++)
461 srenew(inprod[i],snew_size);
463 inprod[noutvec][nframes]=t;
464 /* calculate x: a fitted struture of the selected atoms */
466 reset_x(nfit,ifit,nat,NULL,xread,w_rls);
467 do_fit(nat,w_rls,xref,xread);
469 for (i=0; i<natoms; i++)
470 copy_rvec(xread[index[i]],x[i]);
472 for(v=0; v<noutvec; v++) {
474 /* calculate (mass-weighted) projection */
476 for (i=0; i<natoms; i++) {
477 inp+=(eigvec[vec][i][0]*(x[i][0]-xav[i][0])+
478 eigvec[vec][i][1]*(x[i][1]-xav[i][1])+
479 eigvec[vec][i][2]*(x[i][2]-xav[i][2]))*sqrtm[i];
481 inprod[v][nframes]=inp;
484 for(i=0; i<natoms; i++)
485 for(d=0; d<DIM; d++) {
486 /* misuse xread for output */
487 xread[index[i]][d] = xav[i][d];
488 for(v=0; v<noutvec; v++)
489 xread[index[i]][d] +=
490 inprod[v][nframes]*eigvec[outvec[v]][i][d]/sqrtm[i];
492 write_trx(out,natoms,index,atoms,0,t,box,xread,NULL,NULL);
497 } while (read_next_x(oenv,status,&t,nat,xread,box));
504 snew(xread,atoms->nr);
507 gmx_rmpbc_done(gpbc);
511 snew(ylabel,noutvec);
512 for(v=0; v<noutvec; v++) {
513 sprintf(str,"vec %d",eignr[outvec[v]]+1);
514 ylabel[v]=strdup(str);
516 sprintf(str,"projection on eigenvectors (%s)",proj_unit);
517 write_xvgr_graphs(projfile, noutvec, 1, str, NULL, output_env_get_xvgr_tlabel(oenv),
518 (const char **)ylabel,
519 nframes, inprod[noutvec], inprod, NULL,
520 output_env_get_time_factor(oenv), FALSE, bSplit,oenv);
524 sprintf(str,"projection on eigenvector %d (%s)",
525 eignr[outvec[0]]+1,proj_unit);
526 sprintf(str2,"projection on eigenvector %d (%s)",
527 eignr[outvec[noutvec-1]]+1,proj_unit);
528 xvgrout=xvgropen(twodplotfile,"2D projection of trajectory",str,str2,
530 for(i=0; i<nframes; i++) {
531 if ( bSplit && i>0 && abs(inprod[noutvec][i])<1e-5 )
532 fprintf(xvgrout,"&\n");
533 fprintf(xvgrout,"%10.5f %10.5f\n",inprod[0][i],inprod[noutvec-1][i]);
538 if (threedplotfile) {
543 char *resnm,*atnm, pdbform[STRLEN];
548 gmx_fatal(FARGS,"You have selected less than 3 eigenvectors");
551 bPDB = fn2ftp(threedplotfile)==efPDB;
553 box[XX][XX] = box[YY][YY] = box[ZZ][ZZ] = 1;
555 b4D = bPDB && (noutvec >= 4);
557 fprintf(stderr, "You have selected four or more eigenvectors:\n"
558 "fourth eigenvector will be plotted "
559 "in bfactor field of pdb file\n");
560 sprintf(str,"4D proj. of traj. on eigenv. %d, %d, %d and %d",
561 eignr[outvec[0]]+1,eignr[outvec[1]]+1,
562 eignr[outvec[2]]+1,eignr[outvec[3]]+1);
564 sprintf(str,"3D proj. of traj. on eigenv. %d, %d and %d",
565 eignr[outvec[0]]+1,eignr[outvec[1]]+1,eignr[outvec[2]]+1);
567 init_t_atoms(&atoms,nframes,FALSE);
574 fact=10000.0/nframes;
578 for(i=0; i<nframes; i++) {
579 atoms.atomname[i] = &atnm;
580 atoms.atom[i].resind = i;
581 atoms.resinfo[i].name = &resnm;
582 atoms.resinfo[i].nr = ceil(i*fact);
583 atoms.resinfo[i].ic = ' ';
584 x[i][XX]=inprod[0][i];
585 x[i][YY]=inprod[1][i];
586 x[i][ZZ]=inprod[2][i];
590 if ( ( b4D || bSplit ) && bPDB ) {
591 strcpy(pdbform,pdbformat);
592 strcat(pdbform,"%8.4f%8.4f\n");
594 out=ffopen(threedplotfile,"w");
595 fprintf(out,"HEADER %s\n",str);
597 fprintf(out,"REMARK %s\n","fourth dimension plotted as B-factor");
599 for(i=0; i<atoms.nr; i++) {
600 if ( j>0 && bSplit && abs(inprod[noutvec][i])<1e-5 ) {
601 fprintf(out,"TER\n");
604 fprintf(out,pdbform,"ATOM",i+1,"C","PRJ",' ',j+1,
605 PR_VEC(10*x[i]), 1.0, 10*b[i]);
607 fprintf(out,"CONECT%5d%5d\n", i, i+1);
610 fprintf(out,"TER\n");
613 write_sto_conf(threedplotfile,str,&atoms,x,NULL,ePBC,box);
614 free_t_atoms(&atoms,FALSE);
619 fprintf(stderr,"%11s %17s %17s\n","eigenvector","Minimum","Maximum");
621 "%11s %10s %10s %10s %10s\n","","value","time","value","time");
622 snew(imin,noutvec_extr);
623 snew(imax,noutvec_extr);
624 for(v=0; v<noutvec_extr; v++) {
625 for(i=0; i<nframes; i++) {
626 if (inprod[v][i]<inprod[v][imin[v]])
628 if (inprod[v][i]>inprod[v][imax[v]])
631 min=inprod[v][imin[v]];
632 max=inprod[v][imax[v]];
633 fprintf(stderr,"%7d %10.6f %10.1f %10.6f %10.1f\n",
635 min,inprod[noutvec][imin[v]],max,inprod[noutvec][imax[v]]);
642 /* build format string for filename: */
643 strcpy(str,extremefile);/* copy filename */
644 c=strrchr(str,'.'); /* find where extention begins */
645 strcpy(str2,c); /* get extention */
646 sprintf(c,"%%d%s",str2); /* append '%s' and extention to filename */
647 for(v=0; v<noutvec_extr; v++) {
648 /* make filename using format string */
650 strcpy(str2,extremefile);
652 sprintf(str2,str,eignr[outvec[v]]+1);
653 fprintf(stderr,"Writing %d frames along eigenvector %d to %s\n",
654 nextr,outvec[v]+1,str2);
655 out=open_trx(str2,"w");
656 for(frame=0; frame<nextr; frame++) {
657 if ((extreme==0) && (nextr<=3))
658 for(i=0; i<natoms; i++) {
659 atoms->resinfo[atoms->atom[index[i]].resind].chainid = 'A' + frame;
661 for(i=0; i<natoms; i++)
664 (xav[i][d] + (min*(nextr-frame-1)+max*frame)/(nextr-1)
665 *eigvec[outvec[v]][i][d]/sqrtm[i]);
666 write_trx(out,natoms,index,atoms,0,frame,topbox,xread,NULL,NULL);
671 fprintf(stderr,"\n");
674 static void components(const char *outfile,int natoms,
675 int *eignr,rvec **eigvec,
676 int noutvec,int *outvec,
677 const output_env_t oenv)
681 char str[STRLEN],**ylabel;
683 fprintf(stderr,"Writing eigenvector components to %s\n",outfile);
685 snew(ylabel,noutvec);
688 for(i=0; i<natoms; i++)
690 for(g=0; g<noutvec; g++) {
692 sprintf(str,"vec %d",eignr[v]+1);
693 ylabel[g]=strdup(str);
696 snew(y[g][s],natoms);
698 for(i=0; i<natoms; i++) {
699 y[g][0][i] = norm(eigvec[v][i]);
701 y[g][s+1][i] = eigvec[v][i][s];
704 write_xvgr_graphs(outfile,noutvec,4,"Eigenvector components",
705 "black: total, red: x, green: y, blue: z",
706 "Atom number",(const char **)ylabel,
707 natoms,x,NULL,y,1,FALSE,FALSE,oenv);
708 fprintf(stderr,"\n");
711 static void rmsf(const char *outfile,int natoms,real *sqrtm,
712 int *eignr,rvec **eigvec,
713 int noutvec,int *outvec,
714 real *eigval, int neig,
715 const output_env_t oenv)
719 char str[STRLEN],**ylabel;
721 for(i=0; i<neig; i++)
725 fprintf(stderr,"Writing rmsf to %s\n",outfile);
727 snew(ylabel,noutvec);
730 for(i=0; i<natoms; i++)
732 for(g=0; g<noutvec; g++) {
734 if (eignr[v] >= neig)
735 gmx_fatal(FARGS,"Selected vector %d is larger than the number of eigenvalues (%d)",eignr[v]+1,neig);
736 sprintf(str,"vec %d",eignr[v]+1);
737 ylabel[g]=strdup(str);
739 for(i=0; i<natoms; i++)
740 y[g][i] = sqrt(eigval[eignr[v]]*norm2(eigvec[v][i]))/sqrtm[i];
742 write_xvgr_graphs(outfile,noutvec,1,"RMS fluctuation (nm) ",NULL,
743 "Atom number",(const char **)ylabel,
744 natoms,x,y,NULL,1,TRUE,FALSE,oenv);
745 fprintf(stderr,"\n");
748 int gmx_anaeig(int argc,char *argv[])
750 static const char *desc[] = {
751 "[TT]g_anaeig[tt] analyzes eigenvectors. The eigenvectors can be of a",
752 "covariance matrix ([TT]g_covar[tt]) or of a Normal Modes analysis",
753 "([TT]g_nmeig[tt]).[PAR]",
755 "When a trajectory is projected on eigenvectors, all structures are",
756 "fitted to the structure in the eigenvector file, if present, otherwise",
757 "to the structure in the structure file. When no run input file is",
758 "supplied, periodicity will not be taken into account. Most analyses",
759 "are performed on eigenvectors [TT]-first[tt] to [TT]-last[tt], but when",
760 "[TT]-first[tt] is set to -1 you will be prompted for a selection.[PAR]",
762 "[TT]-comp[tt]: plot the vector components per atom of eigenvectors",
763 "[TT]-first[tt] to [TT]-last[tt].[PAR]",
765 "[TT]-rmsf[tt]: plot the RMS fluctuation per atom of eigenvectors",
766 "[TT]-first[tt] to [TT]-last[tt] (requires [TT]-eig[tt]).[PAR]",
768 "[TT]-proj[tt]: calculate projections of a trajectory on eigenvectors",
769 "[TT]-first[tt] to [TT]-last[tt].",
770 "The projections of a trajectory on the eigenvectors of its",
771 "covariance matrix are called principal components (pc's).",
772 "It is often useful to check the cosine content of the pc's,",
773 "since the pc's of random diffusion are cosines with the number",
774 "of periods equal to half the pc index.",
775 "The cosine content of the pc's can be calculated with the program",
776 "[TT]g_analyze[tt].[PAR]",
778 "[TT]-2d[tt]: calculate a 2d projection of a trajectory on eigenvectors",
779 "[TT]-first[tt] and [TT]-last[tt].[PAR]",
781 "[TT]-3d[tt]: calculate a 3d projection of a trajectory on the first",
782 "three selected eigenvectors.[PAR]",
784 "[TT]-filt[tt]: filter the trajectory to show only the motion along",
785 "eigenvectors [TT]-first[tt] to [TT]-last[tt].[PAR]",
787 "[TT]-extr[tt]: calculate the two extreme projections along a trajectory",
788 "on the average structure and interpolate [TT]-nframes[tt] frames",
789 "between them, or set your own extremes with [TT]-max[tt]. The",
790 "eigenvector [TT]-first[tt] will be written unless [TT]-first[tt] and",
791 "[TT]-last[tt] have been set explicitly, in which case all eigenvectors",
792 "will be written to separate files. Chain identifiers will be added",
793 "when writing a [TT].pdb[tt] file with two or three structures (you",
794 "can use [TT]rasmol -nmrpdb[tt] to view such a [TT].pdb[tt] file).[PAR]",
796 " Overlap calculations between covariance analysis:[BR]",
797 " [BB]Note:[bb] the analysis should use the same fitting structure[PAR]",
799 "[TT]-over[tt]: calculate the subspace overlap of the eigenvectors in",
800 "file [TT]-v2[tt] with eigenvectors [TT]-first[tt] to [TT]-last[tt]",
801 "in file [TT]-v[tt].[PAR]",
803 "[TT]-inpr[tt]: calculate a matrix of inner-products between",
804 "eigenvectors in files [TT]-v[tt] and [TT]-v2[tt]. All eigenvectors",
805 "of both files will be used unless [TT]-first[tt] and [TT]-last[tt]",
806 "have been set explicitly.[PAR]",
808 "When [TT]-v[tt], [TT]-eig[tt], [TT]-v2[tt] and [TT]-eig2[tt] are given,",
809 "a single number for the overlap between the covariance matrices is",
810 "generated. The formulas are:[BR]",
811 " difference = sqrt(tr((sqrt(M1) - sqrt(M2))^2))[BR]",
812 "normalized overlap = 1 - difference/sqrt(tr(M1) + tr(M2))[BR]",
813 " shape overlap = 1 - sqrt(tr((sqrt(M1/tr(M1)) - sqrt(M2/tr(M2)))^2))[BR]",
814 "where M1 and M2 are the two covariance matrices and tr is the trace",
815 "of a matrix. The numbers are proportional to the overlap of the square",
816 "root of the fluctuations. The normalized overlap is the most useful",
817 "number, it is 1 for identical matrices and 0 when the sampled",
818 "subspaces are orthogonal.[PAR]",
819 "When the [TT]-entropy[tt] flag is given an entropy estimate will be",
820 "computed based on the Quasiharmonic approach and based on",
821 "Schlitter's formula."
823 static int first=1,last=8,skip=1,nextr=2,nskip=6;
824 static real max=0.0,temp=298.15;
825 static gmx_bool bSplit=FALSE,bEntropy=FALSE;
827 { "-first", FALSE, etINT, {&first},
828 "First eigenvector for analysis (-1 is select)" },
829 { "-last", FALSE, etINT, {&last},
830 "Last eigenvector for analysis (-1 is till the last)" },
831 { "-skip", FALSE, etINT, {&skip},
832 "Only analyse every nr-th frame" },
833 { "-max", FALSE, etREAL, {&max},
834 "Maximum for projection of the eigenvector on the average structure, "
835 "max=0 gives the extremes" },
836 { "-nframes", FALSE, etINT, {&nextr},
837 "Number of frames for the extremes output" },
838 { "-split", FALSE, etBOOL, {&bSplit},
839 "Split eigenvector projections where time is zero" },
840 { "-entropy", FALSE, etBOOL, {&bEntropy},
841 "Compute entropy according to the Quasiharmonic formula or Schlitter's method." },
842 { "-temp", FALSE, etREAL, {&temp},
843 "Temperature for entropy calculations" },
844 { "-nevskip", FALSE, etINT, {&nskip},
845 "Number of eigenvalues to skip when computing the entropy due to the quasi harmonic approximation. When you do a rotational and/or translational fit prior to the covariance analysis, you get 3 or 6 eigenvalues that are very close to zero, and which should not be taken into account when computing the entropy." }
847 #define NPA asize(pa)
854 rvec *xtop,*xref1,*xref2,*xrefp=NULL;
855 gmx_bool bDMR1,bDMA1,bDMR2,bDMA2;
856 int nvec1,nvec2,*eignr1=NULL,*eignr2=NULL;
857 rvec *x,*xread,*xav1,*xav2,**eigvec1=NULL,**eigvec2=NULL;
859 real xid,totmass,*sqrtm,*w_rls,t,lambda;
862 const char *indexfile;
865 int nout,*iout,noutvec,*outvec,nfit;
866 atom_id *index,*ifit;
867 const char *VecFile,*Vec2File,*topfile;
868 const char *EigFile,*Eig2File;
869 const char *CompFile,*RmsfFile,*ProjOnVecFile;
870 const char *TwoDPlotFile,*ThreeDPlotFile;
871 const char *FilterFile,*ExtremeFile;
872 const char *OverlapFile,*InpMatFile;
873 gmx_bool bFit1,bFit2,bM,bIndex,bTPS,bTop,bVec2,bProj;
874 gmx_bool bFirstToLast,bFirstLastSet,bTraj,bCompare,bPDB3D;
875 real *eigval1=NULL,*eigval2=NULL;
882 { efTRN, "-v", "eigenvec", ffREAD },
883 { efTRN, "-v2", "eigenvec2", ffOPTRD },
884 { efTRX, "-f", NULL, ffOPTRD },
885 { efTPS, NULL, NULL, ffOPTRD },
886 { efNDX, NULL, NULL, ffOPTRD },
887 { efXVG, "-eig", "eigenval", ffOPTRD },
888 { efXVG, "-eig2", "eigenval2", ffOPTRD },
889 { efXVG, "-comp", "eigcomp", ffOPTWR },
890 { efXVG, "-rmsf", "eigrmsf", ffOPTWR },
891 { efXVG, "-proj", "proj", ffOPTWR },
892 { efXVG, "-2d", "2dproj", ffOPTWR },
893 { efSTO, "-3d", "3dproj.pdb", ffOPTWR },
894 { efTRX, "-filt", "filtered", ffOPTWR },
895 { efTRX, "-extr", "extreme.pdb", ffOPTWR },
896 { efXVG, "-over", "overlap", ffOPTWR },
897 { efXPM, "-inpr", "inprod", ffOPTWR }
899 #define NFILE asize(fnm)
901 CopyRight(stderr,argv[0]);
902 parse_common_args(&argc,argv,
903 PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW | PCA_BE_NICE ,
904 NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
906 indexfile=ftp2fn_null(efNDX,NFILE,fnm);
908 VecFile = opt2fn("-v",NFILE,fnm);
909 Vec2File = opt2fn_null("-v2",NFILE,fnm);
910 topfile = ftp2fn(efTPS,NFILE,fnm);
911 EigFile = opt2fn_null("-eig",NFILE,fnm);
912 Eig2File = opt2fn_null("-eig2",NFILE,fnm);
913 CompFile = opt2fn_null("-comp",NFILE,fnm);
914 RmsfFile = opt2fn_null("-rmsf",NFILE,fnm);
915 ProjOnVecFile = opt2fn_null("-proj",NFILE,fnm);
916 TwoDPlotFile = opt2fn_null("-2d",NFILE,fnm);
917 ThreeDPlotFile = opt2fn_null("-3d",NFILE,fnm);
918 FilterFile = opt2fn_null("-filt",NFILE,fnm);
919 ExtremeFile = opt2fn_null("-extr",NFILE,fnm);
920 OverlapFile = opt2fn_null("-over",NFILE,fnm);
921 InpMatFile = ftp2fn_null(efXPM,NFILE,fnm);
923 bTop = fn2bTPX(topfile);
924 bProj = ProjOnVecFile || TwoDPlotFile || ThreeDPlotFile
925 || FilterFile || ExtremeFile;
927 opt2parg_bSet("-first",NPA,pa) && opt2parg_bSet("-last",NPA,pa);
928 bFirstToLast = CompFile || RmsfFile || ProjOnVecFile || FilterFile ||
929 OverlapFile || ((ExtremeFile || InpMatFile) && bFirstLastSet);
930 bVec2 = Vec2File || OverlapFile || InpMatFile;
931 bM = RmsfFile || bProj;
932 bTraj = ProjOnVecFile || FilterFile || (ExtremeFile && (max==0))
933 || TwoDPlotFile || ThreeDPlotFile;
934 bIndex = bM || bProj;
935 bTPS = ftp2bSet(efTPS,NFILE,fnm) || bM || bTraj ||
936 FilterFile || (bIndex && indexfile);
937 bCompare = Vec2File || Eig2File;
938 bPDB3D = fn2ftp(ThreeDPlotFile)==efPDB;
940 read_eigenvectors(VecFile,&natoms,&bFit1,
941 &xref1,&bDMR1,&xav1,&bDMA1,
942 &nvec1,&eignr1,&eigvec1,&eigval1);
945 /* Overwrite eigenvalues from separate files if the user provides them */
946 if (EigFile != NULL) {
947 int neig_tmp = read_xvg(EigFile,&xvgdata,&i);
948 if (neig_tmp != neig1)
949 fprintf(stderr,"Warning: number of eigenvalues in xvg file (%d) does not mtch trr file (%d)\n",neig1,natoms);
951 srenew(eigval1,neig1);
952 for(j=0;j<neig1;j++) {
953 real tmp = eigval1[j];
954 eigval1[j]=xvgdata[1][j];
955 if (debug && (eigval1[j] != tmp))
956 fprintf(debug,"Replacing eigenvalue %d. From trr: %10g, from xvg: %10g\n",
962 fprintf(stderr,"Read %d eigenvalues from %s\n",neig1,EigFile);
966 calc_entropy_qh(stdout,neig1,eigval1,temp,nskip);
967 calc_entropy_schlitter(stdout,neig1,nskip,eigval1,temp);
972 gmx_fatal(FARGS,"Need a second eigenvector file to do this analysis.");
973 read_eigenvectors(Vec2File,&neig2,&bFit2,
974 &xref2,&bDMR2,&xav2,&bDMA2,&nvec2,&eignr2,&eigvec2,&eigval2);
978 gmx_fatal(FARGS,"Dimensions in the eigenvector files don't match");
981 if(Eig2File != NULL) {
982 neig2 = read_xvg(Eig2File,&xvgdata,&i);
983 srenew(eigval2,neig2);
985 eigval2[j]=xvgdata[1][j];
989 fprintf(stderr,"Read %d eigenvalues from %s\n",neig2,Eig2File);
993 if ((!bFit1 || xref1) && !bDMR1 && !bDMA1)
995 if ((xref1==NULL) && (bM || bTraj))
1006 bTop=read_tps_conf(ftp2fn(efTPS,NFILE,fnm),
1007 title,&top,&ePBC,&xtop,NULL,topbox,bM);
1009 gpbc = gmx_rmpbc_init(&top.idef,ePBC,atoms->nr,topbox);
1010 gmx_rmpbc(gpbc,atoms->nr,topbox,xtop);
1011 /* Fitting is only required for the projection */
1012 if (bProj && bFit1) {
1013 if (xref1 == NULL) {
1014 printf("\nNote: the structure in %s should be the same\n"
1015 " as the one used for the fit in g_covar\n",topfile);
1017 printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
1018 get_index(atoms,indexfile,1,&nfit,&ifit,&grpname);
1020 snew(w_rls,atoms->nr);
1021 for(i=0; (i<nfit); i++) {
1023 w_rls[ifit[i]] = atoms->atom[ifit[i]].m;
1025 w_rls[ifit[i]] = 1.0;
1029 snew(xrefp,atoms->nr);
1030 if (xref1 != NULL) {
1031 for(i=0; (i<nfit); i++) {
1032 copy_rvec(xref1[i],xrefp[ifit[i]]);
1035 /* The top coordinates are the fitting reference */
1036 for(i=0; (i<nfit); i++) {
1037 copy_rvec(xtop[ifit[i]],xrefp[ifit[i]]);
1039 reset_x(nfit,ifit,atoms->nr,NULL,xrefp,w_rls);
1042 gmx_rmpbc_done(gpbc);
1046 printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n",natoms);
1047 get_index(atoms,indexfile,1,&i,&index,&grpname);
1049 gmx_fatal(FARGS,"you selected a group with %d elements instead of %d",i,natoms);
1055 proj_unit="u\\S1/2\\Nnm";
1056 for(i=0; (i<natoms); i++)
1057 sqrtm[i]=sqrt(atoms->atom[index[i]].m);
1061 for(i=0; (i<natoms); i++)
1068 for(i=0; (i<natoms); i++)
1069 for(d=0;(d<DIM);d++) {
1070 t+=sqr((xav1[i][d]-xav2[i][d])*sqrtm[i]);
1071 totmass+=sqr(sqrtm[i]);
1073 fprintf(stdout,"RMSD (without fit) between the two average structures:"
1074 " %.3f (nm)\n\n",sqrt(t/totmass));
1081 /* make an index from first to last */
1084 for(i=0; i<nout; i++)
1087 else if (ThreeDPlotFile) {
1088 /* make an index of first+(0,1,2) and last */
1089 nout = bPDB3D ? 4 : 3;
1090 nout = min(last-first+1, nout);
1096 iout[nout-1]=last-1;
1099 /* make an index of first and last */
1107 printf("Select eigenvectors for output, end your selection with 0\n");
1113 srenew(iout,nout+1);
1114 if(1 != scanf("%d",&iout[nout]))
1116 gmx_fatal(FARGS,"Error reading user input");
1120 while (iout[nout]>=0);
1124 /* make an index of the eigenvectors which are present */
1127 for(i=0; i<nout; i++)
1130 while ((j<nvec1) && (eignr1[j]!=iout[i]))
1132 if ((j<nvec1) && (eignr1[j]==iout[i]))
1138 fprintf(stderr,"%d eigenvectors selected for output",noutvec);
1141 fprintf(stderr,":");
1142 for(j=0; j<noutvec; j++)
1143 fprintf(stderr," %d",eignr1[outvec[j]]+1);
1145 fprintf(stderr,"\n");
1148 components(CompFile,natoms,eignr1,eigvec1,noutvec,outvec,oenv);
1151 rmsf(RmsfFile,natoms,sqrtm,eignr1,eigvec1,noutvec,outvec,eigval1,
1155 project(bTraj ? opt2fn("-f",NFILE,fnm) : NULL,
1156 bTop ? &top : NULL,ePBC,topbox,
1157 ProjOnVecFile,TwoDPlotFile,ThreeDPlotFile,FilterFile,skip,
1158 ExtremeFile,bFirstLastSet,max,nextr,atoms,natoms,index,
1159 bFit1,xrefp,nfit,ifit,w_rls,
1160 sqrtm,xav1,eignr1,eigvec1,noutvec,outvec,bSplit,
1164 overlap(OverlapFile,natoms,
1165 eigvec1,nvec2,eignr2,eigvec2,noutvec,outvec,oenv);
1168 inprod_matrix(InpMatFile,natoms,
1169 nvec1,eignr1,eigvec1,nvec2,eignr2,eigvec2,
1170 bFirstLastSet,noutvec,outvec);
1173 compare(natoms,nvec1,eigvec1,nvec2,eigvec2,eigval1,neig1,eigval2,neig2);
1176 if (!CompFile && !bProj && !OverlapFile && !InpMatFile &&
1177 !bCompare && !bEntropy)
1179 fprintf(stderr,"\nIf you want some output,"
1180 " set one (or two or ...) of the output file options\n");
1184 view_all(oenv,NFILE, fnm);