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;
416 real t,inp,**inprod=NULL,min=0,max=0;
417 char str[STRLEN],str2[STRLEN],**ylabel,*c;
419 gmx_rmpbc_t gpbc=NULL;
424 noutvec_extr=noutvec;
430 snew(inprod,noutvec+1);
433 fprintf(stderr,"Writing a filtered trajectory to %s using eigenvectors\n",
435 for(i=0; i<noutvec; i++)
436 fprintf(stderr,"%d ",outvec[i]+1);
437 fprintf(stderr,"\n");
438 out=open_trx(filterfile,"w");
443 nat=read_first_x(oenv,&status,trajfile,&t,&xread,box);
445 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);
449 gpbc = gmx_rmpbc_init(&top->idef,ePBC,nat,box);
454 if (nfr % skip == 0) {
456 gmx_rmpbc(gpbc,nat,box,xread);
457 if (nframes>=snew_size) {
459 for(i=0; i<noutvec+1; i++)
460 srenew(inprod[i],snew_size);
462 inprod[noutvec][nframes]=t;
463 /* calculate x: a fitted struture of the selected atoms */
465 reset_x(nfit,ifit,nat,NULL,xread,w_rls);
466 do_fit(nat,w_rls,xref,xread);
468 for (i=0; i<natoms; i++)
469 copy_rvec(xread[index[i]],x[i]);
471 for(v=0; v<noutvec; v++) {
473 /* calculate (mass-weighted) projection */
475 for (i=0; i<natoms; i++) {
476 inp+=(eigvec[vec][i][0]*(x[i][0]-xav[i][0])+
477 eigvec[vec][i][1]*(x[i][1]-xav[i][1])+
478 eigvec[vec][i][2]*(x[i][2]-xav[i][2]))*sqrtm[i];
480 inprod[v][nframes]=inp;
483 for(i=0; i<natoms; i++)
484 for(d=0; d<DIM; d++) {
485 /* misuse xread for output */
486 xread[index[i]][d] = xav[i][d];
487 for(v=0; v<noutvec; v++)
488 xread[index[i]][d] +=
489 inprod[v][nframes]*eigvec[outvec[v]][i][d]/sqrtm[i];
491 write_trx(out,natoms,index,atoms,0,t,box,xread,NULL,NULL);
496 } while (read_next_x(oenv,status,&t,nat,xread,box));
503 snew(xread,atoms->nr);
506 gmx_rmpbc_done(gpbc);
510 snew(ylabel,noutvec);
511 for(v=0; v<noutvec; v++) {
512 sprintf(str,"vec %d",eignr[outvec[v]]+1);
513 ylabel[v]=strdup(str);
515 sprintf(str,"projection on eigenvectors (%s)",proj_unit);
516 write_xvgr_graphs(projfile, noutvec, 1, str, NULL, output_env_get_xvgr_tlabel(oenv),
517 (const char **)ylabel,
518 nframes, inprod[noutvec], inprod, NULL,
519 output_env_get_time_factor(oenv), FALSE, bSplit,oenv);
523 sprintf(str,"projection on eigenvector %d (%s)",
524 eignr[outvec[0]]+1,proj_unit);
525 sprintf(str2,"projection on eigenvector %d (%s)",
526 eignr[outvec[noutvec-1]]+1,proj_unit);
527 xvgrout=xvgropen(twodplotfile,"2D projection of trajectory",str,str2,
529 for(i=0; i<nframes; i++) {
530 if ( bSplit && i>0 && abs(inprod[noutvec][i])<1e-5 )
531 fprintf(xvgrout,"&\n");
532 fprintf(xvgrout,"%10.5f %10.5f\n",inprod[0][i],inprod[noutvec-1][i]);
537 if (threedplotfile) {
542 char *resnm,*atnm, pdbform[STRLEN];
547 gmx_fatal(FARGS,"You have selected less than 3 eigenvectors");
550 bPDB = fn2ftp(threedplotfile)==efPDB;
552 box[XX][XX] = box[YY][YY] = box[ZZ][ZZ] = 1;
554 b4D = bPDB && (noutvec >= 4);
556 fprintf(stderr, "You have selected four or more eigenvectors:\n"
557 "fourth eigenvector will be plotted "
558 "in bfactor field of pdb file\n");
559 sprintf(str,"4D proj. of traj. on eigenv. %d, %d, %d and %d",
560 eignr[outvec[0]]+1,eignr[outvec[1]]+1,
561 eignr[outvec[2]]+1,eignr[outvec[3]]+1);
563 sprintf(str,"3D proj. of traj. on eigenv. %d, %d and %d",
564 eignr[outvec[0]]+1,eignr[outvec[1]]+1,eignr[outvec[2]]+1);
566 init_t_atoms(&atoms,nframes,FALSE);
573 fact=10000.0/nframes;
577 for(i=0; i<nframes; i++) {
578 atoms.atomname[i] = &atnm;
579 atoms.atom[i].resind = i;
580 atoms.resinfo[i].name = &resnm;
581 atoms.resinfo[i].nr = ceil(i*fact);
582 atoms.resinfo[i].ic = ' ';
583 x[i][XX]=inprod[0][i];
584 x[i][YY]=inprod[1][i];
585 x[i][ZZ]=inprod[2][i];
589 if ( ( b4D || bSplit ) && bPDB ) {
590 strcpy(pdbform,get_pdbformat());
591 strcat(pdbform,"%8.4f%8.4f\n");
593 out=ffopen(threedplotfile,"w");
594 fprintf(out,"HEADER %s\n",str);
596 fprintf(out,"REMARK %s\n","fourth dimension plotted as B-factor");
598 for(i=0; i<atoms.nr; i++) {
599 if ( j>0 && bSplit && abs(inprod[noutvec][i])<1e-5 ) {
600 fprintf(out,"TER\n");
603 fprintf(out,pdbform,"ATOM",i+1,"C","PRJ",' ',j+1,
604 PR_VEC(10*x[i]), 1.0, 10*b[i]);
606 fprintf(out,"CONECT%5d%5d\n", i, i+1);
609 fprintf(out,"TER\n");
612 write_sto_conf(threedplotfile,str,&atoms,x,NULL,ePBC,box);
613 free_t_atoms(&atoms,FALSE);
617 snew(pmin,noutvec_extr);
618 snew(pmax,noutvec_extr);
620 fprintf(stderr,"%11s %17s %17s\n","eigenvector","Minimum","Maximum");
622 "%11s %10s %10s %10s %10s\n","","value","frame","value","frame");
625 for(v=0; v<noutvec_extr; v++) {
626 for(i=0; i<nframes; i++) {
627 if (inprod[v][i]<inprod[v][imin])
629 if (inprod[v][i]>inprod[v][imax])
632 pmin[v] = inprod[v][imin];
633 pmax[v] = inprod[v][imax];
634 fprintf(stderr,"%7d %10.6f %10d %10.6f %10d\n",
636 pmin[v],imin,pmax[v],imax);
643 /* build format string for filename: */
644 strcpy(str,extremefile);/* copy filename */
645 c=strrchr(str,'.'); /* find where extention begins */
646 strcpy(str2,c); /* get extention */
647 sprintf(c,"%%d%s",str2); /* append '%s' and extention to filename */
648 for(v=0; v<noutvec_extr; v++) {
649 /* make filename using format string */
651 strcpy(str2,extremefile);
653 sprintf(str2,str,eignr[outvec[v]]+1);
654 fprintf(stderr,"Writing %d frames along eigenvector %d to %s\n",
655 nextr,outvec[v]+1,str2);
656 out=open_trx(str2,"w");
657 for(frame=0; frame<nextr; frame++) {
658 if ((extreme==0) && (nextr<=3))
659 for(i=0; i<natoms; i++) {
660 atoms->resinfo[atoms->atom[index[i]].resind].chainid = 'A' + frame;
662 for(i=0; i<natoms; i++)
665 (xav[i][d] + (pmin[v]*(nextr-frame-1)+pmax[v]*frame)/(nextr-1)
666 *eigvec[outvec[v]][i][d]/sqrtm[i]);
667 write_trx(out,natoms,index,atoms,0,frame,topbox,xread,NULL,NULL);
674 fprintf(stderr,"\n");
677 static void components(const char *outfile,int natoms,
678 int *eignr,rvec **eigvec,
679 int noutvec,int *outvec,
680 const output_env_t oenv)
684 char str[STRLEN],**ylabel;
686 fprintf(stderr,"Writing eigenvector components to %s\n",outfile);
688 snew(ylabel,noutvec);
691 for(i=0; i<natoms; i++)
693 for(g=0; g<noutvec; g++) {
695 sprintf(str,"vec %d",eignr[v]+1);
696 ylabel[g]=strdup(str);
699 snew(y[g][s],natoms);
701 for(i=0; i<natoms; i++) {
702 y[g][0][i] = norm(eigvec[v][i]);
704 y[g][s+1][i] = eigvec[v][i][s];
707 write_xvgr_graphs(outfile,noutvec,4,"Eigenvector components",
708 "black: total, red: x, green: y, blue: z",
709 "Atom number",(const char **)ylabel,
710 natoms,x,NULL,y,1,FALSE,FALSE,oenv);
711 fprintf(stderr,"\n");
714 static void rmsf(const char *outfile,int natoms,real *sqrtm,
715 int *eignr,rvec **eigvec,
716 int noutvec,int *outvec,
717 real *eigval, int neig,
718 const output_env_t oenv)
722 char str[STRLEN],**ylabel;
724 for(i=0; i<neig; i++)
728 fprintf(stderr,"Writing rmsf to %s\n",outfile);
730 snew(ylabel,noutvec);
733 for(i=0; i<natoms; i++)
735 for(g=0; g<noutvec; g++) {
737 if (eignr[v] >= neig)
738 gmx_fatal(FARGS,"Selected vector %d is larger than the number of eigenvalues (%d)",eignr[v]+1,neig);
739 sprintf(str,"vec %d",eignr[v]+1);
740 ylabel[g]=strdup(str);
742 for(i=0; i<natoms; i++)
743 y[g][i] = sqrt(eigval[eignr[v]]*norm2(eigvec[v][i]))/sqrtm[i];
745 write_xvgr_graphs(outfile,noutvec,1,"RMS fluctuation (nm) ",NULL,
746 "Atom number",(const char **)ylabel,
747 natoms,x,y,NULL,1,TRUE,FALSE,oenv);
748 fprintf(stderr,"\n");
751 int gmx_anaeig(int argc,char *argv[])
753 static const char *desc[] = {
754 "[TT]g_anaeig[tt] analyzes eigenvectors. The eigenvectors can be of a",
755 "covariance matrix ([TT]g_covar[tt]) or of a Normal Modes analysis",
756 "([TT]g_nmeig[tt]).[PAR]",
758 "When a trajectory is projected on eigenvectors, all structures are",
759 "fitted to the structure in the eigenvector file, if present, otherwise",
760 "to the structure in the structure file. When no run input file is",
761 "supplied, periodicity will not be taken into account. Most analyses",
762 "are performed on eigenvectors [TT]-first[tt] to [TT]-last[tt], but when",
763 "[TT]-first[tt] is set to -1 you will be prompted for a selection.[PAR]",
765 "[TT]-comp[tt]: plot the vector components per atom of eigenvectors",
766 "[TT]-first[tt] to [TT]-last[tt].[PAR]",
768 "[TT]-rmsf[tt]: plot the RMS fluctuation per atom of eigenvectors",
769 "[TT]-first[tt] to [TT]-last[tt] (requires [TT]-eig[tt]).[PAR]",
771 "[TT]-proj[tt]: calculate projections of a trajectory on eigenvectors",
772 "[TT]-first[tt] to [TT]-last[tt].",
773 "The projections of a trajectory on the eigenvectors of its",
774 "covariance matrix are called principal components (pc's).",
775 "It is often useful to check the cosine content of the pc's,",
776 "since the pc's of random diffusion are cosines with the number",
777 "of periods equal to half the pc index.",
778 "The cosine content of the pc's can be calculated with the program",
779 "[TT]g_analyze[tt].[PAR]",
781 "[TT]-2d[tt]: calculate a 2d projection of a trajectory on eigenvectors",
782 "[TT]-first[tt] and [TT]-last[tt].[PAR]",
784 "[TT]-3d[tt]: calculate a 3d projection of a trajectory on the first",
785 "three selected eigenvectors.[PAR]",
787 "[TT]-filt[tt]: filter the trajectory to show only the motion along",
788 "eigenvectors [TT]-first[tt] to [TT]-last[tt].[PAR]",
790 "[TT]-extr[tt]: calculate the two extreme projections along a trajectory",
791 "on the average structure and interpolate [TT]-nframes[tt] frames",
792 "between them, or set your own extremes with [TT]-max[tt]. The",
793 "eigenvector [TT]-first[tt] will be written unless [TT]-first[tt] and",
794 "[TT]-last[tt] have been set explicitly, in which case all eigenvectors",
795 "will be written to separate files. Chain identifiers will be added",
796 "when writing a [TT].pdb[tt] file with two or three structures (you",
797 "can use [TT]rasmol -nmrpdb[tt] to view such a [TT].pdb[tt] file).[PAR]",
799 " Overlap calculations between covariance analysis:[BR]",
800 " [BB]Note:[bb] the analysis should use the same fitting structure[PAR]",
802 "[TT]-over[tt]: calculate the subspace overlap of the eigenvectors in",
803 "file [TT]-v2[tt] with eigenvectors [TT]-first[tt] to [TT]-last[tt]",
804 "in file [TT]-v[tt].[PAR]",
806 "[TT]-inpr[tt]: calculate a matrix of inner-products between",
807 "eigenvectors in files [TT]-v[tt] and [TT]-v2[tt]. All eigenvectors",
808 "of both files will be used unless [TT]-first[tt] and [TT]-last[tt]",
809 "have been set explicitly.[PAR]",
811 "When [TT]-v[tt], [TT]-eig[tt], [TT]-v2[tt] and [TT]-eig2[tt] are given,",
812 "a single number for the overlap between the covariance matrices is",
813 "generated. The formulas are:[BR]",
814 " difference = sqrt(tr((sqrt(M1) - sqrt(M2))^2))[BR]",
815 "normalized overlap = 1 - difference/sqrt(tr(M1) + tr(M2))[BR]",
816 " shape overlap = 1 - sqrt(tr((sqrt(M1/tr(M1)) - sqrt(M2/tr(M2)))^2))[BR]",
817 "where M1 and M2 are the two covariance matrices and tr is the trace",
818 "of a matrix. The numbers are proportional to the overlap of the square",
819 "root of the fluctuations. The normalized overlap is the most useful",
820 "number, it is 1 for identical matrices and 0 when the sampled",
821 "subspaces are orthogonal.[PAR]",
822 "When the [TT]-entropy[tt] flag is given an entropy estimate will be",
823 "computed based on the Quasiharmonic approach and based on",
824 "Schlitter's formula."
826 static int first=1,last=-1,skip=1,nextr=2,nskip=6;
827 static real max=0.0,temp=298.15;
828 static gmx_bool bSplit=FALSE,bEntropy=FALSE;
830 { "-first", FALSE, etINT, {&first},
831 "First eigenvector for analysis (-1 is select)" },
832 { "-last", FALSE, etINT, {&last},
833 "Last eigenvector for analysis (-1 is till the last)" },
834 { "-skip", FALSE, etINT, {&skip},
835 "Only analyse every nr-th frame" },
836 { "-max", FALSE, etREAL, {&max},
837 "Maximum for projection of the eigenvector on the average structure, "
838 "max=0 gives the extremes" },
839 { "-nframes", FALSE, etINT, {&nextr},
840 "Number of frames for the extremes output" },
841 { "-split", FALSE, etBOOL, {&bSplit},
842 "Split eigenvector projections where time is zero" },
843 { "-entropy", FALSE, etBOOL, {&bEntropy},
844 "Compute entropy according to the Quasiharmonic formula or Schlitter's method." },
845 { "-temp", FALSE, etREAL, {&temp},
846 "Temperature for entropy calculations" },
847 { "-nevskip", FALSE, etINT, {&nskip},
848 "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." }
850 #define NPA asize(pa)
857 rvec *xtop,*xref1,*xref2,*xrefp=NULL;
858 gmx_bool bDMR1,bDMA1,bDMR2,bDMA2;
859 int nvec1,nvec2,*eignr1=NULL,*eignr2=NULL;
860 rvec *x,*xread,*xav1,*xav2,**eigvec1=NULL,**eigvec2=NULL;
862 real xid,totmass,*sqrtm,*w_rls,t,lambda;
865 const char *indexfile;
868 int nout,*iout,noutvec,*outvec,nfit;
869 atom_id *index,*ifit;
870 const char *VecFile,*Vec2File,*topfile;
871 const char *EigFile,*Eig2File;
872 const char *CompFile,*RmsfFile,*ProjOnVecFile;
873 const char *TwoDPlotFile,*ThreeDPlotFile;
874 const char *FilterFile,*ExtremeFile;
875 const char *OverlapFile,*InpMatFile;
876 gmx_bool bFit1,bFit2,bM,bIndex,bTPS,bTop,bVec2,bProj;
877 gmx_bool bFirstToLast,bFirstLastSet,bTraj,bCompare,bPDB3D;
878 real *eigval1=NULL,*eigval2=NULL;
885 { efTRN, "-v", "eigenvec", ffREAD },
886 { efTRN, "-v2", "eigenvec2", ffOPTRD },
887 { efTRX, "-f", NULL, ffOPTRD },
888 { efTPS, NULL, NULL, ffOPTRD },
889 { efNDX, NULL, NULL, ffOPTRD },
890 { efXVG, "-eig", "eigenval", ffOPTRD },
891 { efXVG, "-eig2", "eigenval2", ffOPTRD },
892 { efXVG, "-comp", "eigcomp", ffOPTWR },
893 { efXVG, "-rmsf", "eigrmsf", ffOPTWR },
894 { efXVG, "-proj", "proj", ffOPTWR },
895 { efXVG, "-2d", "2dproj", ffOPTWR },
896 { efSTO, "-3d", "3dproj.pdb", ffOPTWR },
897 { efTRX, "-filt", "filtered", ffOPTWR },
898 { efTRX, "-extr", "extreme.pdb", ffOPTWR },
899 { efXVG, "-over", "overlap", ffOPTWR },
900 { efXPM, "-inpr", "inprod", ffOPTWR }
902 #define NFILE asize(fnm)
904 CopyRight(stderr,argv[0]);
905 parse_common_args(&argc,argv,
906 PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW | PCA_BE_NICE ,
907 NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
909 indexfile=ftp2fn_null(efNDX,NFILE,fnm);
911 VecFile = opt2fn("-v",NFILE,fnm);
912 Vec2File = opt2fn_null("-v2",NFILE,fnm);
913 topfile = ftp2fn(efTPS,NFILE,fnm);
914 EigFile = opt2fn_null("-eig",NFILE,fnm);
915 Eig2File = opt2fn_null("-eig2",NFILE,fnm);
916 CompFile = opt2fn_null("-comp",NFILE,fnm);
917 RmsfFile = opt2fn_null("-rmsf",NFILE,fnm);
918 ProjOnVecFile = opt2fn_null("-proj",NFILE,fnm);
919 TwoDPlotFile = opt2fn_null("-2d",NFILE,fnm);
920 ThreeDPlotFile = opt2fn_null("-3d",NFILE,fnm);
921 FilterFile = opt2fn_null("-filt",NFILE,fnm);
922 ExtremeFile = opt2fn_null("-extr",NFILE,fnm);
923 OverlapFile = opt2fn_null("-over",NFILE,fnm);
924 InpMatFile = ftp2fn_null(efXPM,NFILE,fnm);
926 bTop = fn2bTPX(topfile);
927 bProj = ProjOnVecFile || TwoDPlotFile || ThreeDPlotFile
928 || FilterFile || ExtremeFile;
930 opt2parg_bSet("-first",NPA,pa) && opt2parg_bSet("-last",NPA,pa);
931 bFirstToLast = CompFile || RmsfFile || ProjOnVecFile || FilterFile ||
932 OverlapFile || ((ExtremeFile || InpMatFile) && bFirstLastSet);
933 bVec2 = Vec2File || OverlapFile || InpMatFile;
934 bM = RmsfFile || bProj;
935 bTraj = ProjOnVecFile || FilterFile || (ExtremeFile && (max==0))
936 || TwoDPlotFile || ThreeDPlotFile;
937 bIndex = bM || bProj;
938 bTPS = ftp2bSet(efTPS,NFILE,fnm) || bM || bTraj ||
939 FilterFile || (bIndex && indexfile);
940 bCompare = Vec2File || Eig2File;
941 bPDB3D = fn2ftp(ThreeDPlotFile)==efPDB;
943 read_eigenvectors(VecFile,&natoms,&bFit1,
944 &xref1,&bDMR1,&xav1,&bDMA1,
945 &nvec1,&eignr1,&eigvec1,&eigval1);
948 /* Overwrite eigenvalues from separate files if the user provides them */
949 if (EigFile != NULL) {
950 int neig_tmp = read_xvg(EigFile,&xvgdata,&i);
951 if (neig_tmp != neig1)
952 fprintf(stderr,"Warning: number of eigenvalues in xvg file (%d) does not mtch trr file (%d)\n",neig1,natoms);
954 srenew(eigval1,neig1);
955 for(j=0;j<neig1;j++) {
956 real tmp = eigval1[j];
957 eigval1[j]=xvgdata[1][j];
958 if (debug && (eigval1[j] != tmp))
959 fprintf(debug,"Replacing eigenvalue %d. From trr: %10g, from xvg: %10g\n",
965 fprintf(stderr,"Read %d eigenvalues from %s\n",neig1,EigFile);
970 gmx_fatal(FARGS,"Can not calculate entropies from mass-weighted eigenvalues, redo the analysis without mass-weighting");
972 calc_entropy_qh(stdout,neig1,eigval1,temp,nskip);
973 calc_entropy_schlitter(stdout,neig1,nskip,eigval1,temp);
978 gmx_fatal(FARGS,"Need a second eigenvector file to do this analysis.");
979 read_eigenvectors(Vec2File,&neig2,&bFit2,
980 &xref2,&bDMR2,&xav2,&bDMA2,&nvec2,&eignr2,&eigvec2,&eigval2);
984 gmx_fatal(FARGS,"Dimensions in the eigenvector files don't match");
987 if(Eig2File != NULL) {
988 neig2 = read_xvg(Eig2File,&xvgdata,&i);
989 srenew(eigval2,neig2);
991 eigval2[j]=xvgdata[1][j];
995 fprintf(stderr,"Read %d eigenvalues from %s\n",neig2,Eig2File);
999 if ((!bFit1 || xref1) && !bDMR1 && !bDMA1)
1001 if ((xref1==NULL) && (bM || bTraj))
1012 bTop=read_tps_conf(ftp2fn(efTPS,NFILE,fnm),
1013 title,&top,&ePBC,&xtop,NULL,topbox,bM);
1015 gpbc = gmx_rmpbc_init(&top.idef,ePBC,atoms->nr,topbox);
1016 gmx_rmpbc(gpbc,atoms->nr,topbox,xtop);
1017 /* Fitting is only required for the projection */
1018 if (bProj && bFit1) {
1019 if (xref1 == NULL) {
1020 printf("\nNote: the structure in %s should be the same\n"
1021 " as the one used for the fit in g_covar\n",topfile);
1023 printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
1024 get_index(atoms,indexfile,1,&nfit,&ifit,&grpname);
1026 snew(w_rls,atoms->nr);
1027 for(i=0; (i<nfit); i++) {
1029 w_rls[ifit[i]] = atoms->atom[ifit[i]].m;
1031 w_rls[ifit[i]] = 1.0;
1035 snew(xrefp,atoms->nr);
1036 if (xref1 != NULL) {
1037 for(i=0; (i<nfit); i++) {
1038 copy_rvec(xref1[i],xrefp[ifit[i]]);
1041 /* The top coordinates are the fitting reference */
1042 for(i=0; (i<nfit); i++) {
1043 copy_rvec(xtop[ifit[i]],xrefp[ifit[i]]);
1045 reset_x(nfit,ifit,atoms->nr,NULL,xrefp,w_rls);
1048 gmx_rmpbc_done(gpbc);
1052 printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n",natoms);
1053 get_index(atoms,indexfile,1,&i,&index,&grpname);
1055 gmx_fatal(FARGS,"you selected a group with %d elements instead of %d",i,natoms);
1061 proj_unit="u\\S1/2\\Nnm";
1062 for(i=0; (i<natoms); i++)
1063 sqrtm[i]=sqrt(atoms->atom[index[i]].m);
1067 for(i=0; (i<natoms); i++)
1074 for(i=0; (i<natoms); i++)
1075 for(d=0;(d<DIM);d++) {
1076 t+=sqr((xav1[i][d]-xav2[i][d])*sqrtm[i]);
1077 totmass+=sqr(sqrtm[i]);
1079 fprintf(stdout,"RMSD (without fit) between the two average structures:"
1080 " %.3f (nm)\n\n",sqrt(t/totmass));
1087 /* make an index from first to last */
1090 for(i=0; i<nout; i++)
1093 else if (ThreeDPlotFile) {
1094 /* make an index of first+(0,1,2) and last */
1095 nout = bPDB3D ? 4 : 3;
1096 nout = min(last-first+1, nout);
1102 iout[nout-1]=last-1;
1105 /* make an index of first and last */
1113 printf("Select eigenvectors for output, end your selection with 0\n");
1119 srenew(iout,nout+1);
1120 if(1 != scanf("%d",&iout[nout]))
1122 gmx_fatal(FARGS,"Error reading user input");
1126 while (iout[nout]>=0);
1130 /* make an index of the eigenvectors which are present */
1133 for(i=0; i<nout; i++)
1136 while ((j<nvec1) && (eignr1[j]!=iout[i]))
1138 if ((j<nvec1) && (eignr1[j]==iout[i]))
1144 fprintf(stderr,"%d eigenvectors selected for output",noutvec);
1147 fprintf(stderr,":");
1148 for(j=0; j<noutvec; j++)
1149 fprintf(stderr," %d",eignr1[outvec[j]]+1);
1151 fprintf(stderr,"\n");
1154 components(CompFile,natoms,eignr1,eigvec1,noutvec,outvec,oenv);
1157 rmsf(RmsfFile,natoms,sqrtm,eignr1,eigvec1,noutvec,outvec,eigval1,
1161 project(bTraj ? opt2fn("-f",NFILE,fnm) : NULL,
1162 bTop ? &top : NULL,ePBC,topbox,
1163 ProjOnVecFile,TwoDPlotFile,ThreeDPlotFile,FilterFile,skip,
1164 ExtremeFile,bFirstLastSet,max,nextr,atoms,natoms,index,
1165 bFit1,xrefp,nfit,ifit,w_rls,
1166 sqrtm,xav1,eignr1,eigvec1,noutvec,outvec,bSplit,
1170 overlap(OverlapFile,natoms,
1171 eigvec1,nvec2,eignr2,eigvec2,noutvec,outvec,oenv);
1174 inprod_matrix(InpMatFile,natoms,
1175 nvec1,eignr1,eigvec1,nvec2,eignr2,eigvec2,
1176 bFirstLastSet,noutvec,outvec);
1179 compare(natoms,nvec1,eigvec1,nvec2,eigvec2,eigval1,neig1,eigval2,neig2);
1182 if (!CompFile && !bProj && !OverlapFile && !InpMatFile &&
1183 !bCompare && !bEntropy)
1185 fprintf(stderr,"\nIf you want some output,"
1186 " set one (or two or ...) of the output file options\n");
1190 view_all(oenv,NFILE, fnm);