2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team,
6 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012,2013, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
48 #include "gmx_fatal.h"
68 static void calc_entropy_qh(FILE *fp,int n,real eigval[],real temp,int nskip)
74 hbar = PLANCK1/(2*M_PI);
75 for(i=0; (i<n-nskip); i++) {
77 lambda = eigval[i]*AMU;
78 w = sqrt(BOLTZMANN*temp/lambda)/NANO;
79 hwkT = (hbar*w)/(BOLTZMANN*temp);
80 dS = (hwkT/(exp(hwkT) - 1) - log(1-exp(-hwkT)));
83 fprintf(debug,"i = %5d w = %10g lam = %10g hwkT = %10g dS = %10g\n",
87 fprintf(stderr,"eigval[%d] = %g\n",i,eigval[i]);
91 fprintf(fp,"The Entropy due to the Quasi Harmonic approximation is %g J/mol K\n",
95 static void calc_entropy_schlitter(FILE *fp,int n,int nskip,
96 real *eigval,real temp)
102 double hbar,kt,kteh,S;
104 hbar = PLANCK1/(2*M_PI);
106 kteh = kt*exp(2.0)/(hbar*hbar)*AMU*(NANO*NANO);
108 fprintf(debug,"n = %d, nskip = %d kteh = %g\n",n,nskip,kteh);
111 for(i=0; (i<n-nskip); i++) {
112 dd = 1+kteh*eigval[i];
117 fprintf(fp,"The Entropy due to the Schlitter formula is %g J/mol K\n",S);
120 const char *proj_unit;
122 static real tick_spacing(real range,int minticks)
130 sp = 0.2*exp(log(10)*ceil(log(range)/log(10)));
131 while (range/sp < minticks-1)
137 static void write_xvgr_graphs(const char *file, int ngraphs, int nsetspergraph,
138 const char *title, const char *subtitle,
139 const char *xlabel, const char **ylabel,
140 int n, real *x, real **y, real ***sy,
141 real scale_x, gmx_bool bZero, gmx_bool bSplit,
142 const output_env_t oenv)
146 real min,max,xsp,ysp;
148 out=ffopen(file,"w");
149 if (output_env_get_xvg_format(oenv) == exvgXMGRACE) {
150 fprintf(out,"@ autoscale onread none\n");
152 for(g=0; g<ngraphs; g++) {
157 if (y[g][i]<min) min=y[g][i];
158 if (y[g][i]>max) max=y[g][i];
163 for(s=0; s<nsetspergraph; s++)
165 if (sy[g][s][i]<min) min=sy[g][s][i];
166 if (sy[g][s][i]>max) max=sy[g][s][i];
172 min=min-0.1*(max-min);
173 max=max+0.1*(max-min);
174 xsp=tick_spacing((x[n-1]-x[0])*scale_x,4);
175 ysp=tick_spacing(max-min,3);
176 if (output_env_get_print_xvgr_codes(oenv)) {
177 fprintf(out,"@ with g%d\n@ g%d on\n",g,g);
179 fprintf(out,"@ title \"%s\"\n",title);
181 fprintf(out,"@ subtitle \"%s\"\n",subtitle);
184 fprintf(out,"@ xaxis label \"%s\"\n",xlabel);
186 fprintf(out,"@ xaxis ticklabel off\n");
188 fprintf(out,"@ world xmin %g\n",x[0]*scale_x);
189 fprintf(out,"@ world xmax %g\n",x[n-1]*scale_x);
190 fprintf(out,"@ world ymin %g\n",min);
191 fprintf(out,"@ world ymax %g\n",max);
193 fprintf(out,"@ view xmin 0.15\n");
194 fprintf(out,"@ view xmax 0.85\n");
195 fprintf(out,"@ view ymin %g\n",0.15+(ngraphs-1-g)*0.7/ngraphs);
196 fprintf(out,"@ view ymax %g\n",0.15+(ngraphs-g)*0.7/ngraphs);
197 fprintf(out,"@ yaxis label \"%s\"\n",ylabel[g]);
198 fprintf(out,"@ xaxis tick major %g\n",xsp);
199 fprintf(out,"@ xaxis tick minor %g\n",xsp/2);
200 fprintf(out,"@ xaxis ticklabel start type spec\n");
201 fprintf(out,"@ xaxis ticklabel start %g\n",ceil(min/xsp)*xsp);
202 fprintf(out,"@ yaxis tick major %g\n",ysp);
203 fprintf(out,"@ yaxis tick minor %g\n",ysp/2);
204 fprintf(out,"@ yaxis ticklabel start type spec\n");
205 fprintf(out,"@ yaxis ticklabel start %g\n",ceil(min/ysp)*ysp);
206 if ((min<0) && (max>0)) {
207 fprintf(out,"@ zeroxaxis bar on\n");
208 fprintf(out,"@ zeroxaxis bar linestyle 3\n");
211 for(s=0; s<nsetspergraph; s++) {
213 if ( bSplit && i>0 && abs(x[i])<1e-5 )
215 if (output_env_get_print_xvgr_codes(oenv))
218 fprintf(out,"%10.4f %10.5f\n",
219 x[i]*scale_x,y ? y[g][i] : sy[g][s][i]);
221 if (output_env_get_print_xvgr_codes(oenv))
229 compare(int natoms,int n1,rvec **eigvec1,int n2,rvec **eigvec2,
230 real *eigval1, int neig1, real *eigval2, int neig2)
234 double sum1,sum2,trace1,trace2,sab,samsb2,tmp,ip;
238 n = min(n,min(neig1,neig2));
239 fprintf(stdout,"Will compare the covariance matrices using %d dimensions\n",n);
247 eigval1[i] = sqrt(eigval1[i]);
250 for(i=n; i<neig1; i++)
251 trace1 += eigval1[i];
258 eigval2[i] = sqrt(eigval2[i]);
261 for(i=n; i<neig2; i++)
262 trace2 += eigval2[i];
264 fprintf(stdout,"Trace of the two matrices: %g and %g\n",sum1,sum2);
265 if (neig1!=n || neig2!=n)
266 fprintf(stdout,"this is %d%% and %d%% of the total trace\n",
267 (int)(100*sum1/trace1+0.5),(int)(100*sum2/trace2+0.5));
268 fprintf(stdout,"Square root of the traces: %g and %g\n",
269 sqrt(sum1),sqrt(sum2));
278 for(k=0; k<natoms; k++)
279 ip += iprod(eigvec1[i][k],eigvec2[j][k]);
280 tmp += eigval2[j]*ip*ip;
282 sab += eigval1[i]*tmp;
285 samsb2 = sum1+sum2-2*sab;
289 fprintf(stdout,"The overlap of the covariance matrices:\n");
290 fprintf(stdout," normalized: %.3f\n",1-sqrt(samsb2/(sum1+sum2)));
291 tmp = 1-sab/sqrt(sum1*sum2);
294 fprintf(stdout," shape: %.3f\n",1-sqrt(tmp));
298 static void inprod_matrix(const char *matfile,int natoms,
299 int nvec1,int *eignr1,rvec **eigvec1,
300 int nvec2,int *eignr2,rvec **eigvec2,
301 gmx_bool bSelect,int noutvec,int *outvec)
305 int i,x1,y1,x,y,nlevels;
307 real inp,*t_x,*t_y,max;
314 for(y1=0; y1<nx; y1++)
315 if (outvec[y1] < nvec2) {
316 t_y[ny] = eignr2[outvec[y1]]+1;
323 t_y[y] = eignr2[y]+1;
326 fprintf(stderr,"Calculating inner-product matrix of %dx%d eigenvectors\n",
332 for(x1=0; x1<nx; x1++) {
338 t_x[x1] = eignr1[x]+1;
339 fprintf(stderr," %d",eignr1[x]+1);
340 for(y1=0; y1<ny; y1++) {
343 while (outvec[y1] >= nvec2)
348 for(i=0; i<natoms; i++)
349 inp += iprod(eigvec1[x][i],eigvec2[y][i]);
350 mat[x1][y1] = fabs(inp);
355 fprintf(stderr,"\n");
356 rlo.r = 1; rlo.g = 1; rlo.b = 1;
357 rhi.r = 0; rhi.g = 0; rhi.b = 0;
359 out = ffopen(matfile,"w");
360 write_xpm(out,0,"Eigenvector inner-products","in.prod.","run 1","run 2",
361 nx,ny,t_x,t_y,mat,0.0,max,rlo,rhi,&nlevels);
365 static void overlap(const char *outfile,int natoms,
367 int nvec2,int *eignr2,rvec **eigvec2,
368 int noutvec,int *outvec,
369 const output_env_t oenv)
375 fprintf(stderr,"Calculating overlap between eigenvectors of set 2 with eigenvectors\n");
376 for(i=0; i<noutvec; i++)
377 fprintf(stderr,"%d ",outvec[i]+1);
378 fprintf(stderr,"\n");
380 out=xvgropen(outfile,"Subspace overlap",
381 "Eigenvectors of trajectory 2","Overlap",oenv);
382 fprintf(out,"@ subtitle \"using %d eigenvectors of trajectory 1\"\n",
385 for(x=0; x<nvec2; x++) {
386 for(v=0; v<noutvec; v++) {
389 for(i=0; i<natoms; i++)
390 inp+=iprod(eigvec1[vec][i],eigvec2[x][i]);
393 fprintf(out,"%5d %5.3f\n",eignr2[x]+1,overlap/noutvec);
399 static void project(const char *trajfile,t_topology *top,int ePBC,matrix topbox,
400 const char *projfile,const char *twodplotfile,
401 const char *threedplotfile, const char *filterfile,int skip,
402 const char *extremefile,gmx_bool bExtrAll,real extreme,
403 int nextr, t_atoms *atoms,int natoms,atom_id *index,
404 gmx_bool bFit,rvec *xref,int nfit,atom_id *ifit,real *w_rls,
405 real *sqrtm,rvec *xav,
406 int *eignr,rvec **eigvec,
407 int noutvec,int *outvec, gmx_bool bSplit,
408 const output_env_t oenv)
411 int nat,i,j,d,v,vec,nfr,nframes=0,snew_size,frame;
412 t_trxstatus *out=NULL;
414 int noutvec_extr,imin,imax;
419 real t,inp,**inprod=NULL,min=0,max=0;
420 char str[STRLEN],str2[STRLEN],**ylabel,*c;
422 gmx_rmpbc_t gpbc=NULL;
427 noutvec_extr=noutvec;
433 snew(inprod,noutvec+1);
436 fprintf(stderr,"Writing a filtered trajectory to %s using eigenvectors\n",
438 for(i=0; i<noutvec; i++)
439 fprintf(stderr,"%d ",outvec[i]+1);
440 fprintf(stderr,"\n");
441 out=open_trx(filterfile,"w");
446 nat=read_first_x(oenv,&status,trajfile,&t,&xread,box);
448 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);
452 gpbc = gmx_rmpbc_init(&top->idef,ePBC,nat,box);
457 if (nfr % skip == 0) {
459 gmx_rmpbc(gpbc,nat,box,xread);
460 if (nframes>=snew_size) {
462 for(i=0; i<noutvec+1; i++)
463 srenew(inprod[i],snew_size);
465 inprod[noutvec][nframes]=t;
466 /* calculate x: a fitted struture of the selected atoms */
468 reset_x(nfit,ifit,nat,NULL,xread,w_rls);
469 do_fit(nat,w_rls,xref,xread);
471 for (i=0; i<natoms; i++)
472 copy_rvec(xread[index[i]],x[i]);
474 for(v=0; v<noutvec; v++) {
476 /* calculate (mass-weighted) projection */
478 for (i=0; i<natoms; i++) {
479 inp+=(eigvec[vec][i][0]*(x[i][0]-xav[i][0])+
480 eigvec[vec][i][1]*(x[i][1]-xav[i][1])+
481 eigvec[vec][i][2]*(x[i][2]-xav[i][2]))*sqrtm[i];
483 inprod[v][nframes]=inp;
486 for(i=0; i<natoms; i++)
487 for(d=0; d<DIM; d++) {
488 /* misuse xread for output */
489 xread[index[i]][d] = xav[i][d];
490 for(v=0; v<noutvec; v++)
491 xread[index[i]][d] +=
492 inprod[v][nframes]*eigvec[outvec[v]][i][d]/sqrtm[i];
494 write_trx(out,natoms,index,atoms,0,t,box,xread,NULL,NULL);
499 } while (read_next_x(oenv,status,&t,nat,xread,box));
506 snew(xread,atoms->nr);
509 gmx_rmpbc_done(gpbc);
513 snew(ylabel,noutvec);
514 for(v=0; v<noutvec; v++) {
515 sprintf(str,"vec %d",eignr[outvec[v]]+1);
516 ylabel[v]=strdup(str);
518 sprintf(str,"projection on eigenvectors (%s)",proj_unit);
519 write_xvgr_graphs(projfile, noutvec, 1, str, NULL, output_env_get_xvgr_tlabel(oenv),
520 (const char **)ylabel,
521 nframes, inprod[noutvec], inprod, NULL,
522 output_env_get_time_factor(oenv), FALSE, bSplit,oenv);
526 sprintf(str,"projection on eigenvector %d (%s)",
527 eignr[outvec[0]]+1,proj_unit);
528 sprintf(str2,"projection on eigenvector %d (%s)",
529 eignr[outvec[noutvec-1]]+1,proj_unit);
530 xvgrout=xvgropen(twodplotfile,"2D projection of trajectory",str,str2,
532 for(i=0; i<nframes; i++) {
533 if ( bSplit && i>0 && abs(inprod[noutvec][i])<1e-5 )
534 fprintf(xvgrout,"&\n");
535 fprintf(xvgrout,"%10.5f %10.5f\n",inprod[0][i],inprod[noutvec-1][i]);
540 if (threedplotfile) {
545 char *resnm,*atnm, pdbform[STRLEN];
550 gmx_fatal(FARGS,"You have selected less than 3 eigenvectors");
553 bPDB = fn2ftp(threedplotfile)==efPDB;
555 box[XX][XX] = box[YY][YY] = box[ZZ][ZZ] = 1;
557 b4D = bPDB && (noutvec >= 4);
559 fprintf(stderr, "You have selected four or more eigenvectors:\n"
560 "fourth eigenvector will be plotted "
561 "in bfactor field of pdb file\n");
562 sprintf(str,"4D proj. of traj. on eigenv. %d, %d, %d and %d",
563 eignr[outvec[0]]+1,eignr[outvec[1]]+1,
564 eignr[outvec[2]]+1,eignr[outvec[3]]+1);
566 sprintf(str,"3D proj. of traj. on eigenv. %d, %d and %d",
567 eignr[outvec[0]]+1,eignr[outvec[1]]+1,eignr[outvec[2]]+1);
569 init_t_atoms(&atoms,nframes,FALSE);
576 fact=10000.0/nframes;
580 for(i=0; i<nframes; i++) {
581 atoms.atomname[i] = &atnm;
582 atoms.atom[i].resind = i;
583 atoms.resinfo[i].name = &resnm;
584 atoms.resinfo[i].nr = ceil(i*fact);
585 atoms.resinfo[i].ic = ' ';
586 x[i][XX]=inprod[0][i];
587 x[i][YY]=inprod[1][i];
588 x[i][ZZ]=inprod[2][i];
592 if ( ( b4D || bSplit ) && bPDB ) {
593 strcpy(pdbform,pdbformat);
594 strcat(pdbform,"%8.4f%8.4f\n");
596 out=ffopen(threedplotfile,"w");
597 fprintf(out,"HEADER %s\n",str);
599 fprintf(out,"REMARK %s\n","fourth dimension plotted as B-factor");
601 for(i=0; i<atoms.nr; i++) {
602 if ( j>0 && bSplit && abs(inprod[noutvec][i])<1e-5 ) {
603 fprintf(out,"TER\n");
606 fprintf(out,pdbform,"ATOM",i+1,"C","PRJ",' ',j+1,
607 PR_VEC(10*x[i]), 1.0, 10*b[i]);
609 fprintf(out,"CONECT%5d%5d\n", i, i+1);
612 fprintf(out,"TER\n");
615 write_sto_conf(threedplotfile,str,&atoms,x,NULL,ePBC,box);
616 free_t_atoms(&atoms,FALSE);
620 snew(pmin,noutvec_extr);
621 snew(pmax,noutvec_extr);
623 fprintf(stderr,"%11s %17s %17s\n","eigenvector","Minimum","Maximum");
625 "%11s %10s %10s %10s %10s\n","","value","frame","value","frame");
628 for(v=0; v<noutvec_extr; v++) {
629 for(i=0; i<nframes; i++) {
630 if (inprod[v][i]<inprod[v][imin])
632 if (inprod[v][i]>inprod[v][imax])
635 pmin[v] = inprod[v][imin];
636 pmax[v] = inprod[v][imax];
637 fprintf(stderr,"%7d %10.6f %10d %10.6f %10d\n",
639 pmin[v],imin,pmax[v],imax);
646 /* build format string for filename: */
647 strcpy(str,extremefile);/* copy filename */
648 c=strrchr(str,'.'); /* find where extention begins */
649 strcpy(str2,c); /* get extention */
650 sprintf(c,"%%d%s",str2); /* append '%s' and extention to filename */
651 for(v=0; v<noutvec_extr; v++) {
652 /* make filename using format string */
654 strcpy(str2,extremefile);
656 sprintf(str2,str,eignr[outvec[v]]+1);
657 fprintf(stderr,"Writing %d frames along eigenvector %d to %s\n",
658 nextr,outvec[v]+1,str2);
659 out=open_trx(str2,"w");
660 for(frame=0; frame<nextr; frame++) {
661 if ((extreme==0) && (nextr<=3))
662 for(i=0; i<natoms; i++) {
663 atoms->resinfo[atoms->atom[index[i]].resind].chainid = 'A' + frame;
665 for(i=0; i<natoms; i++)
668 (xav[i][d] + (pmin[v]*(nextr-frame-1)+pmax[v]*frame)/(nextr-1)
669 *eigvec[outvec[v]][i][d]/sqrtm[i]);
670 write_trx(out,natoms,index,atoms,0,frame,topbox,xread,NULL,NULL);
677 fprintf(stderr,"\n");
680 static void components(const char *outfile,int natoms,
681 int *eignr,rvec **eigvec,
682 int noutvec,int *outvec,
683 const output_env_t oenv)
687 char str[STRLEN],**ylabel;
689 fprintf(stderr,"Writing eigenvector components to %s\n",outfile);
691 snew(ylabel,noutvec);
694 for(i=0; i<natoms; i++)
696 for(g=0; g<noutvec; g++) {
698 sprintf(str,"vec %d",eignr[v]+1);
699 ylabel[g]=strdup(str);
702 snew(y[g][s],natoms);
704 for(i=0; i<natoms; i++) {
705 y[g][0][i] = norm(eigvec[v][i]);
707 y[g][s+1][i] = eigvec[v][i][s];
710 write_xvgr_graphs(outfile,noutvec,4,"Eigenvector components",
711 "black: total, red: x, green: y, blue: z",
712 "Atom number",(const char **)ylabel,
713 natoms,x,NULL,y,1,FALSE,FALSE,oenv);
714 fprintf(stderr,"\n");
717 static void rmsf(const char *outfile,int natoms,real *sqrtm,
718 int *eignr,rvec **eigvec,
719 int noutvec,int *outvec,
720 real *eigval, int neig,
721 const output_env_t oenv)
725 char str[STRLEN],**ylabel;
727 for(i=0; i<neig; i++)
731 fprintf(stderr,"Writing rmsf to %s\n",outfile);
733 snew(ylabel,noutvec);
736 for(i=0; i<natoms; i++)
738 for(g=0; g<noutvec; g++) {
740 if (eignr[v] >= neig)
741 gmx_fatal(FARGS,"Selected vector %d is larger than the number of eigenvalues (%d)",eignr[v]+1,neig);
742 sprintf(str,"vec %d",eignr[v]+1);
743 ylabel[g]=strdup(str);
745 for(i=0; i<natoms; i++)
746 y[g][i] = sqrt(eigval[eignr[v]]*norm2(eigvec[v][i]))/sqrtm[i];
748 write_xvgr_graphs(outfile,noutvec,1,"RMS fluctuation (nm) ",NULL,
749 "Atom number",(const char **)ylabel,
750 natoms,x,y,NULL,1,TRUE,FALSE,oenv);
751 fprintf(stderr,"\n");
754 int gmx_anaeig(int argc,char *argv[])
756 static const char *desc[] = {
757 "[TT]g_anaeig[tt] analyzes eigenvectors. The eigenvectors can be of a",
758 "covariance matrix ([TT]g_covar[tt]) or of a Normal Modes analysis",
759 "([TT]g_nmeig[tt]).[PAR]",
761 "When a trajectory is projected on eigenvectors, all structures are",
762 "fitted to the structure in the eigenvector file, if present, otherwise",
763 "to the structure in the structure file. When no run input file is",
764 "supplied, periodicity will not be taken into account. Most analyses",
765 "are performed on eigenvectors [TT]-first[tt] to [TT]-last[tt], but when",
766 "[TT]-first[tt] is set to -1 you will be prompted for a selection.[PAR]",
768 "[TT]-comp[tt]: plot the vector components per atom of eigenvectors",
769 "[TT]-first[tt] to [TT]-last[tt].[PAR]",
771 "[TT]-rmsf[tt]: plot the RMS fluctuation per atom of eigenvectors",
772 "[TT]-first[tt] to [TT]-last[tt] (requires [TT]-eig[tt]).[PAR]",
774 "[TT]-proj[tt]: calculate projections of a trajectory on eigenvectors",
775 "[TT]-first[tt] to [TT]-last[tt].",
776 "The projections of a trajectory on the eigenvectors of its",
777 "covariance matrix are called principal components (pc's).",
778 "It is often useful to check the cosine content of the pc's,",
779 "since the pc's of random diffusion are cosines with the number",
780 "of periods equal to half the pc index.",
781 "The cosine content of the pc's can be calculated with the program",
782 "[TT]g_analyze[tt].[PAR]",
784 "[TT]-2d[tt]: calculate a 2d projection of a trajectory on eigenvectors",
785 "[TT]-first[tt] and [TT]-last[tt].[PAR]",
787 "[TT]-3d[tt]: calculate a 3d projection of a trajectory on the first",
788 "three selected eigenvectors.[PAR]",
790 "[TT]-filt[tt]: filter the trajectory to show only the motion along",
791 "eigenvectors [TT]-first[tt] to [TT]-last[tt].[PAR]",
793 "[TT]-extr[tt]: calculate the two extreme projections along a trajectory",
794 "on the average structure and interpolate [TT]-nframes[tt] frames",
795 "between them, or set your own extremes with [TT]-max[tt]. The",
796 "eigenvector [TT]-first[tt] will be written unless [TT]-first[tt] and",
797 "[TT]-last[tt] have been set explicitly, in which case all eigenvectors",
798 "will be written to separate files. Chain identifiers will be added",
799 "when writing a [TT].pdb[tt] file with two or three structures (you",
800 "can use [TT]rasmol -nmrpdb[tt] to view such a [TT].pdb[tt] file).[PAR]",
802 " Overlap calculations between covariance analysis:[BR]",
803 " [BB]Note:[bb] the analysis should use the same fitting structure[PAR]",
805 "[TT]-over[tt]: calculate the subspace overlap of the eigenvectors in",
806 "file [TT]-v2[tt] with eigenvectors [TT]-first[tt] to [TT]-last[tt]",
807 "in file [TT]-v[tt].[PAR]",
809 "[TT]-inpr[tt]: calculate a matrix of inner-products between",
810 "eigenvectors in files [TT]-v[tt] and [TT]-v2[tt]. All eigenvectors",
811 "of both files will be used unless [TT]-first[tt] and [TT]-last[tt]",
812 "have been set explicitly.[PAR]",
814 "When [TT]-v[tt], [TT]-eig[tt], [TT]-v2[tt] and [TT]-eig2[tt] are given,",
815 "a single number for the overlap between the covariance matrices is",
816 "generated. The formulas are:[BR]",
817 " difference = sqrt(tr((sqrt(M1) - sqrt(M2))^2))[BR]",
818 "normalized overlap = 1 - difference/sqrt(tr(M1) + tr(M2))[BR]",
819 " shape overlap = 1 - sqrt(tr((sqrt(M1/tr(M1)) - sqrt(M2/tr(M2)))^2))[BR]",
820 "where M1 and M2 are the two covariance matrices and tr is the trace",
821 "of a matrix. The numbers are proportional to the overlap of the square",
822 "root of the fluctuations. The normalized overlap is the most useful",
823 "number, it is 1 for identical matrices and 0 when the sampled",
824 "subspaces are orthogonal.[PAR]",
825 "When the [TT]-entropy[tt] flag is given an entropy estimate will be",
826 "computed based on the Quasiharmonic approach and based on",
827 "Schlitter's formula."
829 static int first=1,last=-1,skip=1,nextr=2,nskip=6;
830 static real max=0.0,temp=298.15;
831 static gmx_bool bSplit=FALSE,bEntropy=FALSE;
833 { "-first", FALSE, etINT, {&first},
834 "First eigenvector for analysis (-1 is select)" },
835 { "-last", FALSE, etINT, {&last},
836 "Last eigenvector for analysis (-1 is till the last)" },
837 { "-skip", FALSE, etINT, {&skip},
838 "Only analyse every nr-th frame" },
839 { "-max", FALSE, etREAL, {&max},
840 "Maximum for projection of the eigenvector on the average structure, "
841 "max=0 gives the extremes" },
842 { "-nframes", FALSE, etINT, {&nextr},
843 "Number of frames for the extremes output" },
844 { "-split", FALSE, etBOOL, {&bSplit},
845 "Split eigenvector projections where time is zero" },
846 { "-entropy", FALSE, etBOOL, {&bEntropy},
847 "Compute entropy according to the Quasiharmonic formula or Schlitter's method." },
848 { "-temp", FALSE, etREAL, {&temp},
849 "Temperature for entropy calculations" },
850 { "-nevskip", FALSE, etINT, {&nskip},
851 "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." }
853 #define NPA asize(pa)
860 rvec *xtop,*xref1,*xref2,*xrefp=NULL;
861 gmx_bool bDMR1,bDMA1,bDMR2,bDMA2;
862 int nvec1,nvec2,*eignr1=NULL,*eignr2=NULL;
863 rvec *x,*xread,*xav1,*xav2,**eigvec1=NULL,**eigvec2=NULL;
865 real xid,totmass,*sqrtm,*w_rls,t,lambda;
868 const char *indexfile;
871 int nout,*iout,noutvec,*outvec,nfit;
872 atom_id *index,*ifit;
873 const char *VecFile,*Vec2File,*topfile;
874 const char *EigFile,*Eig2File;
875 const char *CompFile,*RmsfFile,*ProjOnVecFile;
876 const char *TwoDPlotFile,*ThreeDPlotFile;
877 const char *FilterFile,*ExtremeFile;
878 const char *OverlapFile,*InpMatFile;
879 gmx_bool bFit1,bFit2,bM,bIndex,bTPS,bTop,bVec2,bProj;
880 gmx_bool bFirstToLast,bFirstLastSet,bTraj,bCompare,bPDB3D;
881 real *eigval1=NULL,*eigval2=NULL;
888 { efTRN, "-v", "eigenvec", ffREAD },
889 { efTRN, "-v2", "eigenvec2", ffOPTRD },
890 { efTRX, "-f", NULL, ffOPTRD },
891 { efTPS, NULL, NULL, ffOPTRD },
892 { efNDX, NULL, NULL, ffOPTRD },
893 { efXVG, "-eig", "eigenval", ffOPTRD },
894 { efXVG, "-eig2", "eigenval2", ffOPTRD },
895 { efXVG, "-comp", "eigcomp", ffOPTWR },
896 { efXVG, "-rmsf", "eigrmsf", ffOPTWR },
897 { efXVG, "-proj", "proj", ffOPTWR },
898 { efXVG, "-2d", "2dproj", ffOPTWR },
899 { efSTO, "-3d", "3dproj.pdb", ffOPTWR },
900 { efTRX, "-filt", "filtered", ffOPTWR },
901 { efTRX, "-extr", "extreme.pdb", ffOPTWR },
902 { efXVG, "-over", "overlap", ffOPTWR },
903 { efXPM, "-inpr", "inprod", ffOPTWR }
905 #define NFILE asize(fnm)
907 CopyRight(stderr,argv[0]);
908 parse_common_args(&argc,argv,
909 PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW | PCA_BE_NICE ,
910 NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
912 indexfile=ftp2fn_null(efNDX,NFILE,fnm);
914 VecFile = opt2fn("-v",NFILE,fnm);
915 Vec2File = opt2fn_null("-v2",NFILE,fnm);
916 topfile = ftp2fn(efTPS,NFILE,fnm);
917 EigFile = opt2fn_null("-eig",NFILE,fnm);
918 Eig2File = opt2fn_null("-eig2",NFILE,fnm);
919 CompFile = opt2fn_null("-comp",NFILE,fnm);
920 RmsfFile = opt2fn_null("-rmsf",NFILE,fnm);
921 ProjOnVecFile = opt2fn_null("-proj",NFILE,fnm);
922 TwoDPlotFile = opt2fn_null("-2d",NFILE,fnm);
923 ThreeDPlotFile = opt2fn_null("-3d",NFILE,fnm);
924 FilterFile = opt2fn_null("-filt",NFILE,fnm);
925 ExtremeFile = opt2fn_null("-extr",NFILE,fnm);
926 OverlapFile = opt2fn_null("-over",NFILE,fnm);
927 InpMatFile = ftp2fn_null(efXPM,NFILE,fnm);
929 bTop = fn2bTPX(topfile);
930 bProj = ProjOnVecFile || TwoDPlotFile || ThreeDPlotFile
931 || FilterFile || ExtremeFile;
933 opt2parg_bSet("-first",NPA,pa) && opt2parg_bSet("-last",NPA,pa);
934 bFirstToLast = CompFile || RmsfFile || ProjOnVecFile || FilterFile ||
935 OverlapFile || ((ExtremeFile || InpMatFile) && bFirstLastSet);
936 bVec2 = Vec2File || OverlapFile || InpMatFile;
937 bM = RmsfFile || bProj;
938 bTraj = ProjOnVecFile || FilterFile || (ExtremeFile && (max==0))
939 || TwoDPlotFile || ThreeDPlotFile;
940 bIndex = bM || bProj;
941 bTPS = ftp2bSet(efTPS,NFILE,fnm) || bM || bTraj ||
942 FilterFile || (bIndex && indexfile);
943 bCompare = Vec2File || Eig2File;
944 bPDB3D = fn2ftp(ThreeDPlotFile)==efPDB;
946 read_eigenvectors(VecFile,&natoms,&bFit1,
947 &xref1,&bDMR1,&xav1,&bDMA1,
948 &nvec1,&eignr1,&eigvec1,&eigval1);
951 /* Overwrite eigenvalues from separate files if the user provides them */
952 if (EigFile != NULL) {
953 int neig_tmp = read_xvg(EigFile,&xvgdata,&i);
954 if (neig_tmp != neig1)
955 fprintf(stderr,"Warning: number of eigenvalues in xvg file (%d) does not mtch trr file (%d)\n",neig1,natoms);
957 srenew(eigval1,neig1);
958 for(j=0;j<neig1;j++) {
959 real tmp = eigval1[j];
960 eigval1[j]=xvgdata[1][j];
961 if (debug && (eigval1[j] != tmp))
962 fprintf(debug,"Replacing eigenvalue %d. From trr: %10g, from xvg: %10g\n",
968 fprintf(stderr,"Read %d eigenvalues from %s\n",neig1,EigFile);
973 gmx_fatal(FARGS,"Can not calculate entropies from mass-weighted eigenvalues, redo the analysis without mass-weighting");
975 calc_entropy_qh(stdout,neig1,eigval1,temp,nskip);
976 calc_entropy_schlitter(stdout,neig1,nskip,eigval1,temp);
981 gmx_fatal(FARGS,"Need a second eigenvector file to do this analysis.");
982 read_eigenvectors(Vec2File,&neig2,&bFit2,
983 &xref2,&bDMR2,&xav2,&bDMA2,&nvec2,&eignr2,&eigvec2,&eigval2);
987 gmx_fatal(FARGS,"Dimensions in the eigenvector files don't match");
990 if(Eig2File != NULL) {
991 neig2 = read_xvg(Eig2File,&xvgdata,&i);
992 srenew(eigval2,neig2);
994 eigval2[j]=xvgdata[1][j];
998 fprintf(stderr,"Read %d eigenvalues from %s\n",neig2,Eig2File);
1002 if ((!bFit1 || xref1) && !bDMR1 && !bDMA1)
1004 if ((xref1==NULL) && (bM || bTraj))
1015 bTop=read_tps_conf(ftp2fn(efTPS,NFILE,fnm),
1016 title,&top,&ePBC,&xtop,NULL,topbox,bM);
1018 gpbc = gmx_rmpbc_init(&top.idef,ePBC,atoms->nr,topbox);
1019 gmx_rmpbc(gpbc,atoms->nr,topbox,xtop);
1020 /* Fitting is only required for the projection */
1021 if (bProj && bFit1) {
1022 if (xref1 == NULL) {
1023 printf("\nNote: the structure in %s should be the same\n"
1024 " as the one used for the fit in g_covar\n",topfile);
1026 printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
1027 get_index(atoms,indexfile,1,&nfit,&ifit,&grpname);
1029 snew(w_rls,atoms->nr);
1030 for(i=0; (i<nfit); i++) {
1032 w_rls[ifit[i]] = atoms->atom[ifit[i]].m;
1034 w_rls[ifit[i]] = 1.0;
1038 snew(xrefp,atoms->nr);
1039 if (xref1 != NULL) {
1040 for(i=0; (i<nfit); i++) {
1041 copy_rvec(xref1[i],xrefp[ifit[i]]);
1044 /* The top coordinates are the fitting reference */
1045 for(i=0; (i<nfit); i++) {
1046 copy_rvec(xtop[ifit[i]],xrefp[ifit[i]]);
1048 reset_x(nfit,ifit,atoms->nr,NULL,xrefp,w_rls);
1051 gmx_rmpbc_done(gpbc);
1055 printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n",natoms);
1056 get_index(atoms,indexfile,1,&i,&index,&grpname);
1058 gmx_fatal(FARGS,"you selected a group with %d elements instead of %d",i,natoms);
1064 proj_unit="u\\S1/2\\Nnm";
1065 for(i=0; (i<natoms); i++)
1066 sqrtm[i]=sqrt(atoms->atom[index[i]].m);
1070 for(i=0; (i<natoms); i++)
1077 for(i=0; (i<natoms); i++)
1078 for(d=0;(d<DIM);d++) {
1079 t+=sqr((xav1[i][d]-xav2[i][d])*sqrtm[i]);
1080 totmass+=sqr(sqrtm[i]);
1082 fprintf(stdout,"RMSD (without fit) between the two average structures:"
1083 " %.3f (nm)\n\n",sqrt(t/totmass));
1090 /* make an index from first to last */
1093 for(i=0; i<nout; i++)
1096 else if (ThreeDPlotFile) {
1097 /* make an index of first+(0,1,2) and last */
1098 nout = bPDB3D ? 4 : 3;
1099 nout = min(last-first+1, nout);
1105 iout[nout-1]=last-1;
1108 /* make an index of first and last */
1116 printf("Select eigenvectors for output, end your selection with 0\n");
1122 srenew(iout,nout+1);
1123 if(1 != scanf("%d",&iout[nout]))
1125 gmx_fatal(FARGS,"Error reading user input");
1129 while (iout[nout]>=0);
1133 /* make an index of the eigenvectors which are present */
1136 for(i=0; i<nout; i++)
1139 while ((j<nvec1) && (eignr1[j]!=iout[i]))
1141 if ((j<nvec1) && (eignr1[j]==iout[i]))
1147 fprintf(stderr,"%d eigenvectors selected for output",noutvec);
1150 fprintf(stderr,":");
1151 for(j=0; j<noutvec; j++)
1152 fprintf(stderr," %d",eignr1[outvec[j]]+1);
1154 fprintf(stderr,"\n");
1157 components(CompFile,natoms,eignr1,eigvec1,noutvec,outvec,oenv);
1160 rmsf(RmsfFile,natoms,sqrtm,eignr1,eigvec1,noutvec,outvec,eigval1,
1164 project(bTraj ? opt2fn("-f",NFILE,fnm) : NULL,
1165 bTop ? &top : NULL,ePBC,topbox,
1166 ProjOnVecFile,TwoDPlotFile,ThreeDPlotFile,FilterFile,skip,
1167 ExtremeFile,bFirstLastSet,max,nextr,atoms,natoms,index,
1168 bFit1,xrefp,nfit,ifit,w_rls,
1169 sqrtm,xav1,eignr1,eigvec1,noutvec,outvec,bSplit,
1173 overlap(OverlapFile,natoms,
1174 eigvec1,nvec2,eignr2,eigvec2,noutvec,outvec,oenv);
1177 inprod_matrix(InpMatFile,natoms,
1178 nvec1,eignr1,eigvec1,nvec2,eignr2,eigvec2,
1179 bFirstLastSet,noutvec,outvec);
1182 compare(natoms,nvec1,eigvec1,nvec2,eigvec2,eigval1,neig1,eigval2,neig2);
1185 if (!CompFile && !bProj && !OverlapFile && !InpMatFile &&
1186 !bCompare && !bEntropy)
1188 fprintf(stderr,"\nIf you want some output,"
1189 " set one (or two or ...) of the output file options\n");
1193 view_all(oenv,NFILE, fnm);