src/tools/gmx_anaeig.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38 #include <math.h>
  39 #include <string.h>
  40 #include "statutil.h"
  41 #include "sysstuff.h"
  42 #include "typedefs.h"
  43 #include "smalloc.h"
  44 #include "macros.h"
  45 #include "gmx_fatal.h"
  46 #include "vec.h"
  47 #include "pbc.h"
  48 #include "copyrite.h"
  49 #include "futil.h"
  50 #include "statutil.h"
  51 #include "index.h"
  52 #include "pdbio.h"
  53 #include "confio.h"
  54 #include "tpxio.h"
  55 #include "matio.h"
  56 #include "mshift.h"
  57 #include "xvgr.h"
  58 #include "do_fit.h"
  59 #include "rmpbc.h"
  60 #include "txtdump.h"
  61 #include "eigio.h"
  62 #include "physics.h"
  63 #include "gmx_ana.h"
  64
  65 static void calc_entropy_qh(FILE *fp,int n,real eigval[],real temp,int nskip)
  66 {
  67   int    i;
  68   double hwkT,w,dS,S=0;
  69   double hbar,lambda;
  70
  71   hbar = PLANCK1/(2*M_PI);
  72   for(i=0; (i<n-nskip); i++) {
  73     if (eigval[i] > 0) {
  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)));
  78       S     += dS;
  79       if (debug)
  80         fprintf(debug,"i = %5d w = %10g lam = %10g hwkT = %10g dS = %10g\n",
  81                 i,w,lambda,hwkT,dS);
  82     }
  83     else {
  84       fprintf(stderr,"eigval[%d] = %g\n",i,eigval[i]);
  85       w = 0;
  86     }
  87   }
  88   fprintf(fp,"The Entropy due to the Quasi Harmonic approximation is %g J/mol K\n",
  89           S*RGAS);
  90 }
  91
  92 static void calc_entropy_schlitter(FILE *fp,int n,int nskip,
  93                                    real *eigval,real temp)
  94 {
  95   double dd,deter;
  96   int    *indx;
  97   int    i,j,k,m;
  98   char   buf[256];
  99   double hbar,kt,kteh,S;
 100
 101   hbar = PLANCK1/(2*M_PI);
 102   kt   = BOLTZMANN*temp;
 103   kteh = kt*exp(2.0)/(hbar*hbar)*AMU*(NANO*NANO);
 104   if (debug)
 105     fprintf(debug,"n = %d, nskip = %d kteh = %g\n",n,nskip,kteh);
 106
 107   deter = 0;
 108   for(i=0; (i<n-nskip); i++) {
 109     dd = 1+kteh*eigval[i];
 110     deter += log(dd);
 111   }
 112   S = 0.5*RGAS*deter;
 113
 114   fprintf(fp,"The Entropy due to the Schlitter formula is %g J/mol K\n",S);
 115 }
 116
 117 const char *proj_unit;
 118
 119 static real tick_spacing(real range,int minticks)
 120 {
 121   real sp;
 122
 123   if (range <= 0) {
 124     return 1.0;
 125   }
 126
 127   sp = 0.2*exp(log(10)*ceil(log(range)/log(10)));
 128   while (range/sp < minticks-1)
 129     sp = sp/2;
 130
 131   return sp;
 132 }
 133
 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)
 140 {
 141   FILE *out;
 142   int g,s,i;
 143   real min,max,xsp,ysp;
 144
 145   out=ffopen(file,"w");
 146   if (output_env_get_xvg_format(oenv) == exvgXMGRACE) {
 147     fprintf(out,"@ autoscale onread none\n");
 148   }
 149   for(g=0; g<ngraphs; g++) {
 150     if (y) {
 151       min=y[g][0];
 152       max=y[g][0];
 153       for(i=0; i<n; i++) {
 154         if (y[g][i]<min) min=y[g][i];
 155         if (y[g][i]>max) max=y[g][i];
 156       }
 157     } else {
 158       min=sy[g][0][0];
 159       max=sy[g][0][0];
 160       for(s=0; s<nsetspergraph; s++)
 161         for(i=0; i<n; i++) {
 162           if (sy[g][s][i]<min) min=sy[g][s][i];
 163           if (sy[g][s][i]>max) max=sy[g][s][i];
 164         }
 165     }
 166     if (bZero)
 167       min=0;
 168     else
 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);
 175       if (g==0) {
 176         fprintf(out,"@ title \"%s\"\n",title);
 177         if (subtitle)
 178           fprintf(out,"@ subtitle \"%s\"\n",subtitle);
 179       }
 180       if (g==ngraphs-1)
 181         fprintf(out,"@ xaxis  label \"%s\"\n",xlabel);
 182       else
 183         fprintf(out,"@ xaxis  ticklabel off\n");
 184       if (n > 1) {
 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);
 189       }
 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");
 206       }
 207     }
 208     for(s=0; s<nsetspergraph; s++) {
 209         for(i=0; i<n; i++) {
 210             if ( bSplit && i>0 && abs(x[i])<1e-5 )
 211             {
 212                 if (output_env_get_print_xvgr_codes(oenv))
 213                     fprintf(out,"&\n");
 214             }
 215             fprintf(out,"%10.4f %10.5f\n",
 216                     x[i]*scale_x,y ? y[g][i] : sy[g][s][i]);
 217         }
 218         if (output_env_get_print_xvgr_codes(oenv))
 219             fprintf(out,"&\n");
 220     }
 221   }
 222   ffclose(out);
 223 }
 224
 225 static void
 226 compare(int natoms,int n1,rvec **eigvec1,int n2,rvec **eigvec2,
 227         real *eigval1, int neig1, real *eigval2, int neig2)
 228 {
 229   int    n,nrow;
 230   int    i,j,k;
 231   double sum1,sum2,trace1,trace2,sab,samsb2,tmp,ip;
 232
 233   n = min(n1,n2);
 234
 235   n = min(n,min(neig1,neig2));
 236   fprintf(stdout,"Will compare the covariance matrices using %d dimensions\n",n);
 237
 238   sum1 = 0;
 239   for(i=0; i<n; i++)
 240   {
 241       if (eigval1[i] < 0)
 242           eigval1[i] = 0;
 243       sum1 += eigval1[i];
 244       eigval1[i] = sqrt(eigval1[i]);
 245   }
 246   trace1 = sum1;
 247   for(i=n; i<neig1; i++)
 248       trace1 += eigval1[i];
 249   sum2 = 0;
 250   for(i=0; i<n; i++)
 251   {
 252     if (eigval2[i] < 0)
 253         eigval2[i] = 0;
 254     sum2 += eigval2[i];
 255     eigval2[i] = sqrt(eigval2[i]);
 256   }
 257   trace2 = sum2;
 258   for(i=n; i<neig2; i++)
 259     trace2 += eigval2[i];
 260
 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));
 267
 268   sab = 0;
 269   for(i=0; i<n; i++)
 270   {
 271     tmp = 0;
 272     for(j=0; j<n; j++)
 273     {
 274       ip = 0;
 275       for(k=0; k<natoms; k++)
 276           ip += iprod(eigvec1[i][k],eigvec2[j][k]);
 277       tmp += eigval2[j]*ip*ip;
 278     }
 279     sab += eigval1[i]*tmp;
 280   }
 281
 282   samsb2 = sum1+sum2-2*sab;
 283   if (samsb2 < 0)
 284       samsb2 = 0;
 285
 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);
 289   if (tmp < 0)
 290       tmp = 0;
 291   fprintf(stdout,"       shape:  %.3f\n",1-sqrt(tmp));
 292 }
 293
 294
 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)
 299 {
 300   FILE  *out;
 301   real  **mat;
 302   int   i,x1,y1,x,y,nlevels;
 303   int   nx,ny;
 304   real  inp,*t_x,*t_y,max;
 305   t_rgb rlo,rhi;
 306
 307   snew(t_y,nvec2);
 308   if (bSelect) {
 309     nx = noutvec;
 310     ny = 0;
 311     for(y1=0; y1<nx; y1++)
 312       if (outvec[y1] < nvec2) {
 313         t_y[ny] = eignr2[outvec[y1]]+1;
 314         ny++;
 315       }
 316   } else {
 317     nx = nvec1;
 318     ny = nvec2;
 319     for(y=0; y<ny; y++)
 320       t_y[y] = eignr2[y]+1;
 321   }
 322
 323   fprintf(stderr,"Calculating inner-product matrix of %dx%d eigenvectors\n",
 324           nx,nvec2);
 325
 326   snew(mat,nx);
 327   snew(t_x,nx);
 328   max = 0;
 329   for(x1=0; x1<nx; x1++) {
 330     snew(mat[x1],ny);
 331     if (bSelect)
 332       x = outvec[x1];
 333     else
 334       x = x1;
 335     t_x[x1] = eignr1[x]+1;
 336     fprintf(stderr," %d",eignr1[x]+1);
 337     for(y1=0; y1<ny; y1++) {
 338       inp = 0;
 339       if (bSelect) {
 340         while (outvec[y1] >= nvec2)
 341           y1++;
 342         y= outvec[y1];
 343       } else
 344         y = y1;
 345       for(i=0; i<natoms; i++)
 346         inp += iprod(eigvec1[x][i],eigvec2[y][i]);
 347       mat[x1][y1] = fabs(inp);
 348       if (mat[x1][y1]>max)
 349         max = mat[x1][y1];
 350     }
 351   }
 352   fprintf(stderr,"\n");
 353   rlo.r = 1; rlo.g = 1; rlo.b = 1;
 354   rhi.r = 0; rhi.g = 0; rhi.b = 0;
 355   nlevels = 41;
 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);
 359   ffclose(out);
 360 }
 361
 362 static void overlap(const char *outfile,int natoms,
 363                     rvec **eigvec1,
 364                     int nvec2,int *eignr2,rvec **eigvec2,
 365                     int noutvec,int *outvec,
 366                     const output_env_t oenv)
 367 {
 368   FILE *out;
 369   int i,v,vec,x;
 370   real overlap,inp;
 371
 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");
 376
 377   out=xvgropen(outfile,"Subspace overlap",
 378                "Eigenvectors of trajectory 2","Overlap",oenv);
 379   fprintf(out,"@ subtitle \"using %d eigenvectors of trajectory 1\"\n",
 380           noutvec);
 381   overlap=0;
 382   for(x=0; x<nvec2; x++) {
 383     for(v=0; v<noutvec; v++) {
 384       vec=outvec[v];
 385       inp=0;
 386       for(i=0; i<natoms; i++)
 387         inp+=iprod(eigvec1[vec][i],eigvec2[x][i]);
 388       overlap+=sqr(inp);
 389     }
 390     fprintf(out,"%5d  %5.3f\n",eignr2[x]+1,overlap/noutvec);
 391   }
 392
 393   ffclose(out);
 394 }
 395
 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)
 406 {
 407   FILE    *xvgrout=NULL;
 408   int     nat,i,j,d,v,vec,nfr,nframes=0,snew_size,frame;
 409   t_trxstatus *out=NULL;
 410   t_trxstatus *status;
 411   int     noutvec_extr,imin,imax;
 412   real    *pmin,*pmax;
 413   atom_id *all_at;
 414   matrix  box;
 415   rvec    *xread,*x;
 416   real    t,inp,**inprod=NULL,min=0,max=0;
 417   char    str[STRLEN],str2[STRLEN],**ylabel,*c;
 418   real    fact;
 419   gmx_rmpbc_t  gpbc=NULL;
 420
 421   snew(x,natoms);
 422
 423   if (bExtrAll)
 424     noutvec_extr=noutvec;
 425   else
 426     noutvec_extr=1;
 427
 428
 429   if (trajfile) {
 430     snew(inprod,noutvec+1);
 431
 432     if (filterfile) {
 433       fprintf(stderr,"Writing a filtered trajectory to %s using eigenvectors\n",
 434               filterfile);
 435       for(i=0; i<noutvec; i++)
 436         fprintf(stderr,"%d ",outvec[i]+1);
 437       fprintf(stderr,"\n");
 438       out=open_trx(filterfile,"w");
 439     }
 440     snew_size=0;
 441     nfr=0;
 442     nframes=0;
 443     nat=read_first_x(oenv,&status,trajfile,&t,&xread,box);
 444     if (nat>atoms->nr)
 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);
 446     snew(all_at,nat);
 447
 448     if (top)
 449       gpbc = gmx_rmpbc_init(&top->idef,ePBC,nat,box);
 450     if (top)
 451     gmx_rmpbc_done(gpbc);
 452
 453     for(i=0; i<nat; i++)
 454       all_at[i]=i;
 455     do {
 456       if (nfr % skip == 0) {
 457         if (top)
 458           gmx_rmpbc(gpbc,nat,box,xread);
 459         if (nframes>=snew_size) {
 460           snew_size+=100;
 461           for(i=0; i<noutvec+1; i++)
 462             srenew(inprod[i],snew_size);
 463         }
 464         inprod[noutvec][nframes]=t;
 465         /* calculate x: a fitted struture of the selected atoms */
 466         if (bFit) {
 467           reset_x(nfit,ifit,nat,NULL,xread,w_rls);
 468           do_fit(nat,w_rls,xref,xread);
 469         }
 470         for (i=0; i<natoms; i++)
 471           copy_rvec(xread[index[i]],x[i]);
 472
 473         for(v=0; v<noutvec; v++) {
 474           vec=outvec[v];
 475           /* calculate (mass-weighted) projection */
 476           inp=0;
 477           for (i=0; i<natoms; i++) {
 478             inp+=(eigvec[vec][i][0]*(x[i][0]-xav[i][0])+
 479             eigvec[vec][i][1]*(x[i][1]-xav[i][1])+
 480             eigvec[vec][i][2]*(x[i][2]-xav[i][2]))*sqrtm[i];
 481           }
 482           inprod[v][nframes]=inp;
 483         }
 484         if (filterfile) {
 485           for(i=0; i<natoms; i++)
 486             for(d=0; d<DIM; d++) {
 487               /* misuse xread for output */
 488               xread[index[i]][d] = xav[i][d];
 489               for(v=0; v<noutvec; v++)
 490                 xread[index[i]][d] +=
 491                   inprod[v][nframes]*eigvec[outvec[v]][i][d]/sqrtm[i];
 492             }
 493           write_trx(out,natoms,index,atoms,0,t,box,xread,NULL,NULL);
 494         }
 495         nframes++;
 496       }
 497       nfr++;
 498     } while (read_next_x(oenv,status,&t,nat,xread,box));
 499     close_trx(status);
 500      sfree(x);
 501      if (filterfile)
 502        close_trx(out);
 503   }
 504   else
 505     snew(xread,atoms->nr);
 506
 507   if (top)
 508     gmx_rmpbc_done(gpbc);
 509
 510
 511   if (projfile) {
 512     snew(ylabel,noutvec);
 513     for(v=0; v<noutvec; v++) {
 514       sprintf(str,"vec %d",eignr[outvec[v]]+1);
 515       ylabel[v]=strdup(str);
 516     }
 517     sprintf(str,"projection on eigenvectors (%s)",proj_unit);
 518     write_xvgr_graphs(projfile, noutvec, 1, str, NULL, output_env_get_xvgr_tlabel(oenv),
 519                       (const char **)ylabel,
 520                       nframes, inprod[noutvec], inprod, NULL,
 521                       output_env_get_time_factor(oenv), FALSE, bSplit,oenv);
 522   }
 523
 524   if (twodplotfile) {
 525     sprintf(str,"projection on eigenvector %d (%s)",
 526             eignr[outvec[0]]+1,proj_unit);
 527     sprintf(str2,"projection on eigenvector %d (%s)",
 528             eignr[outvec[noutvec-1]]+1,proj_unit);
 529     xvgrout=xvgropen(twodplotfile,"2D projection of trajectory",str,str2,
 530                      oenv);
 531     for(i=0; i<nframes; i++) {
 532       if ( bSplit && i>0 && abs(inprod[noutvec][i])<1e-5 )
 533         fprintf(xvgrout,"&\n");
 534       fprintf(xvgrout,"%10.5f %10.5f\n",inprod[0][i],inprod[noutvec-1][i]);
 535     }
 536     ffclose(xvgrout);
 537   }
 538
 539   if (threedplotfile) {
 540     t_atoms atoms;
 541     rvec    *x;
 542     real    *b=NULL;
 543     matrix  box;
 544     char    *resnm,*atnm, pdbform[STRLEN];
 545     gmx_bool    bPDB, b4D;
 546     FILE    *out;
 547
 548     if (noutvec < 3)
 549       gmx_fatal(FARGS,"You have selected less than 3 eigenvectors");
 550
 551     /* initialize */
 552     bPDB = fn2ftp(threedplotfile)==efPDB;
 553     clear_mat(box);
 554     box[XX][XX] = box[YY][YY] = box[ZZ][ZZ] = 1;
 555
 556     b4D = bPDB && (noutvec >= 4);
 557     if (b4D) {
 558       fprintf(stderr, "You have selected four or more eigenvectors:\n"
 559               "fourth eigenvector will be plotted "
 560               "in bfactor field of pdb file\n");
 561       sprintf(str,"4D proj. of traj. on eigenv. %d, %d, %d and %d",
 562               eignr[outvec[0]]+1,eignr[outvec[1]]+1,
 563               eignr[outvec[2]]+1,eignr[outvec[3]]+1);
 564     } else {
 565     sprintf(str,"3D proj. of traj. on eigenv. %d, %d and %d",
 566             eignr[outvec[0]]+1,eignr[outvec[1]]+1,eignr[outvec[2]]+1);
 567     }
 568     init_t_atoms(&atoms,nframes,FALSE);
 569     snew(x,nframes);
 570     snew(b,nframes);
 571     atnm=strdup("C");
 572     resnm=strdup("PRJ");
 573
 574     if(nframes>10000)
 575       fact=10000.0/nframes;
 576     else
 577       fact=1.0;
 578
 579     for(i=0; i<nframes; i++) {
 580       atoms.atomname[i] = &atnm;
 581       atoms.atom[i].resind = i;
 582       atoms.resinfo[i].name = &resnm;
 583       atoms.resinfo[i].nr   = ceil(i*fact);
 584       atoms.resinfo[i].ic   = ' ';
 585       x[i][XX]=inprod[0][i];
 586       x[i][YY]=inprod[1][i];
 587       x[i][ZZ]=inprod[2][i];
 588       if (b4D)
 589         b[i]  =inprod[3][i];
 590     }
 591     if ( ( b4D || bSplit ) && bPDB ) {
 592       strcpy(pdbform,pdbformat);
 593       strcat(pdbform,"%8.4f%8.4f\n");
 594
 595       out=ffopen(threedplotfile,"w");
 596       fprintf(out,"HEADER    %s\n",str);
 597       if ( b4D )
 598         fprintf(out,"REMARK    %s\n","fourth dimension plotted as B-factor");
 599       j=0;
 600       for(i=0; i<atoms.nr; i++) {
 601         if ( j>0 && bSplit && abs(inprod[noutvec][i])<1e-5 ) {
 602           fprintf(out,"TER\n");
 603           j=0;
 604         }
 605         fprintf(out,pdbform,"ATOM",i+1,"C","PRJ",' ',j+1,
 606                 PR_VEC(10*x[i]), 1.0, 10*b[i]);
 607         if (j>0)
 608           fprintf(out,"CONECT%5d%5d\n", i, i+1);
 609         j++;
 610       }
 611       fprintf(out,"TER\n");
 612       ffclose(out);
 613     } else
 614       write_sto_conf(threedplotfile,str,&atoms,x,NULL,ePBC,box);
 615     free_t_atoms(&atoms,FALSE);
 616   }
 617
 618   if (extremefile) {
 619     snew(pmin,noutvec_extr);
 620     snew(pmax,noutvec_extr);
 621     if (extreme==0) {
 622       fprintf(stderr,"%11s %17s %17s\n","eigenvector","Minimum","Maximum");
 623       fprintf(stderr,
 624               "%11s %10s %10s %10s %10s\n","","value","frame","value","frame");
 625       imin = 0;
 626       imax = 0;
 627       for(v=0; v<noutvec_extr; v++) {
 628         for(i=0; i<nframes; i++) {
 629           if (inprod[v][i]<inprod[v][imin])
 630             imin = i;
 631           if (inprod[v][i]>inprod[v][imax])
 632             imax = i;
 633         }
 634         pmin[v] = inprod[v][imin];
 635         pmax[v] = inprod[v][imax];
 636         fprintf(stderr,"%7d     %10.6f %10d %10.6f %10d\n",
 637                 eignr[outvec[v]]+1,
 638                 pmin[v],imin,pmax[v],imax);
 639       }
 640     }
 641     else {
 642       pmin[0] = -extreme;
 643       pmax[0] =  extreme;
 644     }
 645     /* build format string for filename: */
 646     strcpy(str,extremefile);/* copy filename */
 647     c=strrchr(str,'.'); /* find where extention begins */
 648     strcpy(str2,c); /* get extention */
 649     sprintf(c,"%%d%s",str2); /* append '%s' and extention to filename */
 650     for(v=0; v<noutvec_extr; v++) {
 651       /* make filename using format string */
 652       if (noutvec_extr==1)
 653         strcpy(str2,extremefile);
 654       else
 655         sprintf(str2,str,eignr[outvec[v]]+1);
 656       fprintf(stderr,"Writing %d frames along eigenvector %d to %s\n",
 657               nextr,outvec[v]+1,str2);
 658       out=open_trx(str2,"w");
 659       for(frame=0; frame<nextr; frame++) {
 660         if ((extreme==0) && (nextr<=3))
 661           for(i=0; i<natoms; i++) {
 662             atoms->resinfo[atoms->atom[index[i]].resind].chainid = 'A' + frame;
 663           }
 664         for(i=0; i<natoms; i++)
 665           for(d=0; d<DIM; d++)
 666             xread[index[i]][d] =
 667               (xav[i][d] + (pmin[v]*(nextr-frame-1)+pmax[v]*frame)/(nextr-1)
 668               *eigvec[outvec[v]][i][d]/sqrtm[i]);
 669         write_trx(out,natoms,index,atoms,0,frame,topbox,xread,NULL,NULL);
 670       }
 671       close_trx(out);
 672     }
 673     sfree(pmin);
 674     sfree(pmax);
 675   }
 676   fprintf(stderr,"\n");
 677 }
 678
 679 static void components(const char *outfile,int natoms,
 680                        int *eignr,rvec **eigvec,
 681                        int noutvec,int *outvec,
 682                        const output_env_t oenv)
 683 {
 684   int g,s,v,i;
 685   real *x,***y;
 686   char str[STRLEN],**ylabel;
 687
 688   fprintf(stderr,"Writing eigenvector components to %s\n",outfile);
 689
 690   snew(ylabel,noutvec);
 691   snew(y,noutvec);
 692   snew(x,natoms);
 693   for(i=0; i<natoms; i++)
 694     x[i]=i+1;
 695   for(g=0; g<noutvec; g++) {
 696     v=outvec[g];
 697     sprintf(str,"vec %d",eignr[v]+1);
 698     ylabel[g]=strdup(str);
 699     snew(y[g],4);
 700     for(s=0; s<4; s++) {
 701       snew(y[g][s],natoms);
 702     }
 703     for(i=0; i<natoms; i++) {
 704       y[g][0][i] = norm(eigvec[v][i]);
 705       for(s=0; s<3; s++)
 706         y[g][s+1][i] = eigvec[v][i][s];
 707     }
 708   }
 709   write_xvgr_graphs(outfile,noutvec,4,"Eigenvector components",
 710                     "black: total, red: x, green: y, blue: z",
 711                     "Atom number",(const char **)ylabel,
 712                     natoms,x,NULL,y,1,FALSE,FALSE,oenv);
 713   fprintf(stderr,"\n");
 714 }
 715
 716 static void rmsf(const char *outfile,int natoms,real *sqrtm,
 717                  int *eignr,rvec **eigvec,
 718                  int noutvec,int *outvec,
 719                  real *eigval, int neig,
 720                  const output_env_t oenv)
 721 {
 722   int   nrow,g,v,i;
 723   real  *x,**y;
 724   char  str[STRLEN],**ylabel;
 725
 726   for(i=0; i<neig; i++)
 727       if (eigval[i] < 0)
 728           eigval[i] = 0;
 729
 730   fprintf(stderr,"Writing rmsf to %s\n",outfile);
 731
 732   snew(ylabel,noutvec);
 733   snew(y,noutvec);
 734   snew(x,natoms);
 735   for(i=0; i<natoms; i++)
 736       x[i]=i+1;
 737   for(g=0; g<noutvec; g++) {
 738       v=outvec[g];
 739       if (eignr[v] >= neig)
 740           gmx_fatal(FARGS,"Selected vector %d is larger than the number of eigenvalues (%d)",eignr[v]+1,neig);
 741       sprintf(str,"vec %d",eignr[v]+1);
 742       ylabel[g]=strdup(str);
 743       snew(y[g],natoms);
 744       for(i=0; i<natoms; i++)
 745           y[g][i] = sqrt(eigval[eignr[v]]*norm2(eigvec[v][i]))/sqrtm[i];
 746   }
 747   write_xvgr_graphs(outfile,noutvec,1,"RMS fluctuation (nm) ",NULL,
 748                     "Atom number",(const char **)ylabel,
 749                     natoms,x,y,NULL,1,TRUE,FALSE,oenv);
 750   fprintf(stderr,"\n");
 751 }
 752
 753 int gmx_anaeig(int argc,char *argv[])
 754 {
 755   static const char *desc[] = {
 756     "[TT]g_anaeig[tt] analyzes eigenvectors. The eigenvectors can be of a",
 757     "covariance matrix ([TT]g_covar[tt]) or of a Normal Modes analysis",
 758     "([TT]g_nmeig[tt]).[PAR]",
 759
 760     "When a trajectory is projected on eigenvectors, all structures are",
 761     "fitted to the structure in the eigenvector file, if present, otherwise",
 762     "to the structure in the structure file. When no run input file is",
 763     "supplied, periodicity will not be taken into account. Most analyses",
 764     "are performed on eigenvectors [TT]-first[tt] to [TT]-last[tt], but when",
 765     "[TT]-first[tt] is set to -1 you will be prompted for a selection.[PAR]",
 766
 767     "[TT]-comp[tt]: plot the vector components per atom of eigenvectors",
 768     "[TT]-first[tt] to [TT]-last[tt].[PAR]",
 769
 770     "[TT]-rmsf[tt]: plot the RMS fluctuation per atom of eigenvectors",
 771     "[TT]-first[tt] to [TT]-last[tt] (requires [TT]-eig[tt]).[PAR]",
 772
 773     "[TT]-proj[tt]: calculate projections of a trajectory on eigenvectors",
 774     "[TT]-first[tt] to [TT]-last[tt].",
 775     "The projections of a trajectory on the eigenvectors of its",
 776     "covariance matrix are called principal components (pc's).",
 777     "It is often useful to check the cosine content of the pc's,",
 778     "since the pc's of random diffusion are cosines with the number",
 779     "of periods equal to half the pc index.",
 780     "The cosine content of the pc's can be calculated with the program",
 781     "[TT]g_analyze[tt].[PAR]",
 782
 783     "[TT]-2d[tt]: calculate a 2d projection of a trajectory on eigenvectors",
 784     "[TT]-first[tt] and [TT]-last[tt].[PAR]",
 785
 786     "[TT]-3d[tt]: calculate a 3d projection of a trajectory on the first",
 787     "three selected eigenvectors.[PAR]",
 788
 789     "[TT]-filt[tt]: filter the trajectory to show only the motion along",
 790     "eigenvectors [TT]-first[tt] to [TT]-last[tt].[PAR]",
 791
 792     "[TT]-extr[tt]: calculate the two extreme projections along a trajectory",
 793     "on the average structure and interpolate [TT]-nframes[tt] frames",
 794     "between them, or set your own extremes with [TT]-max[tt]. The",
 795     "eigenvector [TT]-first[tt] will be written unless [TT]-first[tt] and",
 796     "[TT]-last[tt] have been set explicitly, in which case all eigenvectors",
 797     "will be written to separate files. Chain identifiers will be added",
 798     "when writing a [TT].pdb[tt] file with two or three structures (you",
 799     "can use [TT]rasmol -nmrpdb[tt] to view such a [TT].pdb[tt] file).[PAR]",
 800
 801     "  Overlap calculations between covariance analysis:[BR]",
 802     "  [BB]Note:[bb] the analysis should use the same fitting structure[PAR]",
 803
 804     "[TT]-over[tt]: calculate the subspace overlap of the eigenvectors in",
 805     "file [TT]-v2[tt] with eigenvectors [TT]-first[tt] to [TT]-last[tt]",
 806     "in file [TT]-v[tt].[PAR]",
 807
 808     "[TT]-inpr[tt]: calculate a matrix of inner-products between",
 809     "eigenvectors in files [TT]-v[tt] and [TT]-v2[tt]. All eigenvectors",
 810     "of both files will be used unless [TT]-first[tt] and [TT]-last[tt]",
 811     "have been set explicitly.[PAR]",
 812
 813     "When [TT]-v[tt], [TT]-eig[tt], [TT]-v2[tt] and [TT]-eig2[tt] are given,",
 814     "a single number for the overlap between the covariance matrices is",
 815     "generated. The formulas are:[BR]",
 816     "        difference = sqrt(tr((sqrt(M1) - sqrt(M2))^2))[BR]",
 817     "normalized overlap = 1 - difference/sqrt(tr(M1) + tr(M2))[BR]",
 818     "     shape overlap = 1 - sqrt(tr((sqrt(M1/tr(M1)) - sqrt(M2/tr(M2)))^2))[BR]",
 819     "where M1 and M2 are the two covariance matrices and tr is the trace",
 820     "of a matrix. The numbers are proportional to the overlap of the square",
 821     "root of the fluctuations. The normalized overlap is the most useful",
 822     "number, it is 1 for identical matrices and 0 when the sampled",
 823     "subspaces are orthogonal.[PAR]",
 824     "When the [TT]-entropy[tt] flag is given an entropy estimate will be",
 825     "computed based on the Quasiharmonic approach and based on",
 826     "Schlitter's formula."
 827   };
 828   static int  first=1,last=-1,skip=1,nextr=2,nskip=6;
 829   static real max=0.0,temp=298.15;
 830   static gmx_bool bSplit=FALSE,bEntropy=FALSE;
 831   t_pargs pa[] = {
 832     { "-first", FALSE, etINT, {&first},
 833       "First eigenvector for analysis (-1 is select)" },
 834     { "-last",  FALSE, etINT, {&last},
 835       "Last eigenvector for analysis (-1 is till the last)" },
 836     { "-skip",  FALSE, etINT, {&skip},
 837       "Only analyse every nr-th frame" },
 838     { "-max",  FALSE, etREAL, {&max},
 839       "Maximum for projection of the eigenvector on the average structure, "
 840       "max=0 gives the extremes" },
 841     { "-nframes",  FALSE, etINT, {&nextr},
 842       "Number of frames for the extremes output" },
 843     { "-split",   FALSE, etBOOL, {&bSplit},
 844       "Split eigenvector projections where time is zero" },
 845     { "-entropy", FALSE, etBOOL, {&bEntropy},
 846       "Compute entropy according to the Quasiharmonic formula or Schlitter's method." },
 847     { "-temp",    FALSE, etREAL, {&temp},
 848       "Temperature for entropy calculations" },
 849     { "-nevskip", FALSE, etINT, {&nskip},
 850       "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." }
 851   };
 852 #define NPA asize(pa)
 853
 854   FILE       *out;
 855   int        status,trjout;
 856   t_topology top;
 857   int        ePBC=-1;
 858   t_atoms    *atoms=NULL;
 859   rvec       *xtop,*xref1,*xref2,*xrefp=NULL;
 860   gmx_bool       bDMR1,bDMA1,bDMR2,bDMA2;
 861   int        nvec1,nvec2,*eignr1=NULL,*eignr2=NULL;
 862   rvec       *x,*xread,*xav1,*xav2,**eigvec1=NULL,**eigvec2=NULL;
 863   matrix     topbox;
 864   real       xid,totmass,*sqrtm,*w_rls,t,lambda;
 865   int        natoms,step;
 866   char       *grpname;
 867   const char *indexfile;
 868   char       title[STRLEN];
 869   int        i,j,d;
 870   int        nout,*iout,noutvec,*outvec,nfit;
 871   atom_id    *index,*ifit;
 872   const char *VecFile,*Vec2File,*topfile;
 873   const char *EigFile,*Eig2File;
 874   const char *CompFile,*RmsfFile,*ProjOnVecFile;
 875   const char *TwoDPlotFile,*ThreeDPlotFile;
 876   const char *FilterFile,*ExtremeFile;
 877   const char *OverlapFile,*InpMatFile;
 878   gmx_bool       bFit1,bFit2,bM,bIndex,bTPS,bTop,bVec2,bProj;
 879   gmx_bool       bFirstToLast,bFirstLastSet,bTraj,bCompare,bPDB3D;
 880   real       *eigval1=NULL,*eigval2=NULL;
 881   int        neig1,neig2;
 882   double     **xvgdata;
 883   output_env_t oenv;
 884   gmx_rmpbc_t  gpbc;
 885
 886   t_filenm fnm[] = {
 887     { efTRN, "-v",    "eigenvec",    ffREAD  },
 888     { efTRN, "-v2",   "eigenvec2",   ffOPTRD },
 889     { efTRX, "-f",    NULL,          ffOPTRD },
 890     { efTPS, NULL,    NULL,          ffOPTRD },
 891     { efNDX, NULL,    NULL,          ffOPTRD },
 892     { efXVG, "-eig", "eigenval",     ffOPTRD },
 893     { efXVG, "-eig2", "eigenval2",   ffOPTRD },
 894     { efXVG, "-comp", "eigcomp",     ffOPTWR },
 895     { efXVG, "-rmsf", "eigrmsf",     ffOPTWR },
 896     { efXVG, "-proj", "proj",        ffOPTWR },
 897     { efXVG, "-2d",   "2dproj",      ffOPTWR },
 898     { efSTO, "-3d",   "3dproj.pdb",  ffOPTWR },
 899     { efTRX, "-filt", "filtered",    ffOPTWR },
 900     { efTRX, "-extr", "extreme.pdb", ffOPTWR },
 901     { efXVG, "-over", "overlap",     ffOPTWR },
 902     { efXPM, "-inpr", "inprod",      ffOPTWR }
 903   };
 904 #define NFILE asize(fnm)
 905
 906   CopyRight(stderr,argv[0]);
 907   parse_common_args(&argc,argv,
 908                     PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW | PCA_BE_NICE ,
 909                     NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
 910
 911   indexfile=ftp2fn_null(efNDX,NFILE,fnm);
 912
 913   VecFile         = opt2fn("-v",NFILE,fnm);
 914   Vec2File        = opt2fn_null("-v2",NFILE,fnm);
 915   topfile         = ftp2fn(efTPS,NFILE,fnm);
 916   EigFile         = opt2fn_null("-eig",NFILE,fnm);
 917   Eig2File        = opt2fn_null("-eig2",NFILE,fnm);
 918   CompFile        = opt2fn_null("-comp",NFILE,fnm);
 919   RmsfFile        = opt2fn_null("-rmsf",NFILE,fnm);
 920   ProjOnVecFile   = opt2fn_null("-proj",NFILE,fnm);
 921   TwoDPlotFile    = opt2fn_null("-2d",NFILE,fnm);
 922   ThreeDPlotFile  = opt2fn_null("-3d",NFILE,fnm);
 923   FilterFile      = opt2fn_null("-filt",NFILE,fnm);
 924   ExtremeFile     = opt2fn_null("-extr",NFILE,fnm);
 925   OverlapFile     = opt2fn_null("-over",NFILE,fnm);
 926   InpMatFile      = ftp2fn_null(efXPM,NFILE,fnm);
 927
 928   bTop   = fn2bTPX(topfile);
 929   bProj  = ProjOnVecFile || TwoDPlotFile || ThreeDPlotFile
 930     || FilterFile || ExtremeFile;
 931   bFirstLastSet  =
 932     opt2parg_bSet("-first",NPA,pa) && opt2parg_bSet("-last",NPA,pa);
 933   bFirstToLast = CompFile || RmsfFile || ProjOnVecFile || FilterFile ||
 934     OverlapFile || ((ExtremeFile || InpMatFile) && bFirstLastSet);
 935   bVec2  = Vec2File || OverlapFile || InpMatFile;
 936   bM     = RmsfFile || bProj;
 937   bTraj  = ProjOnVecFile || FilterFile || (ExtremeFile && (max==0))
 938     || TwoDPlotFile || ThreeDPlotFile;
 939   bIndex = bM || bProj;
 940   bTPS   = ftp2bSet(efTPS,NFILE,fnm) || bM || bTraj ||
 941     FilterFile  || (bIndex && indexfile);
 942   bCompare = Vec2File || Eig2File;
 943   bPDB3D = fn2ftp(ThreeDPlotFile)==efPDB;
 944
 945   read_eigenvectors(VecFile,&natoms,&bFit1,
 946                     &xref1,&bDMR1,&xav1,&bDMA1,
 947                     &nvec1,&eignr1,&eigvec1,&eigval1);
 948   neig1 = DIM*natoms;
 949
 950   /* Overwrite eigenvalues from separate files if the user provides them */
 951   if (EigFile != NULL) {
 952     int neig_tmp = read_xvg(EigFile,&xvgdata,&i);
 953     if (neig_tmp != neig1)
 954       fprintf(stderr,"Warning: number of eigenvalues in xvg file (%d) does not mtch trr file (%d)\n",neig1,natoms);
 955     neig1 = neig_tmp;
 956     srenew(eigval1,neig1);
 957     for(j=0;j<neig1;j++) {
 958       real tmp = eigval1[j];
 959       eigval1[j]=xvgdata[1][j];
 960       if (debug && (eigval1[j] != tmp))
 961         fprintf(debug,"Replacing eigenvalue %d. From trr: %10g, from xvg: %10g\n",
 962                 j,tmp,eigval1[j]);
 963     }
 964     for(j=0;j<i;j++)
 965       sfree(xvgdata[j]);
 966     sfree(xvgdata);
 967     fprintf(stderr,"Read %d eigenvalues from %s\n",neig1,EigFile);
 968   }
 969
 970   if (bEntropy) {
 971     if (bDMA1) {
 972       gmx_fatal(FARGS,"Can not calculate entropies from mass-weighted eigenvalues, redo the analysis without mass-weighting");
 973     }
 974     calc_entropy_qh(stdout,neig1,eigval1,temp,nskip);
 975     calc_entropy_schlitter(stdout,neig1,nskip,eigval1,temp);
 976   }
 977
 978   if (bVec2) {
 979     if (!Vec2File)
 980       gmx_fatal(FARGS,"Need a second eigenvector file to do this analysis.");
 981     read_eigenvectors(Vec2File,&neig2,&bFit2,
 982                       &xref2,&bDMR2,&xav2,&bDMA2,&nvec2,&eignr2,&eigvec2,&eigval2);
 983
 984     neig2 = DIM*neig2;
 985     if (neig2 != neig1)
 986       gmx_fatal(FARGS,"Dimensions in the eigenvector files don't match");
 987   }
 988
 989   if(Eig2File != NULL) {
 990     neig2 = read_xvg(Eig2File,&xvgdata,&i);
 991     srenew(eigval2,neig2);
 992     for(j=0;j<neig2;j++)
 993       eigval2[j]=xvgdata[1][j];
 994     for(j=0;j<i;j++)
 995       sfree(xvgdata[j]);
 996     sfree(xvgdata);
 997     fprintf(stderr,"Read %d eigenvalues from %s\n",neig2,Eig2File);
 998   }
 999
1000
1001   if ((!bFit1 || xref1) && !bDMR1 && !bDMA1)
1002     bM=FALSE;
1003   if ((xref1==NULL) && (bM || bTraj))
1004     bTPS=TRUE;
1005
1006   xtop=NULL;
1007   nfit=0;
1008   ifit=NULL;
1009   w_rls=NULL;
1010
1011   if (!bTPS) {
1012     bTop=FALSE;
1013   } else {
1014     bTop=read_tps_conf(ftp2fn(efTPS,NFILE,fnm),
1015                        title,&top,&ePBC,&xtop,NULL,topbox,bM);
1016     atoms=&top.atoms;
1017     gpbc = gmx_rmpbc_init(&top.idef,ePBC,atoms->nr,topbox);
1018     gmx_rmpbc(gpbc,atoms->nr,topbox,xtop);
1019     /* Fitting is only required for the projection */
1020     if (bProj && bFit1) {
1021       if (xref1 == NULL) {
1022           printf("\nNote: the structure in %s should be the same\n"
1023                  "      as the one used for the fit in g_covar\n",topfile);
1024       }
1025       printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
1026       get_index(atoms,indexfile,1,&nfit,&ifit,&grpname);
1027
1028       snew(w_rls,atoms->nr);
1029       for(i=0; (i<nfit); i++) {
1030         if (bDMR1) {
1031           w_rls[ifit[i]] = atoms->atom[ifit[i]].m;
1032         } else {
1033           w_rls[ifit[i]] = 1.0;
1034         }
1035       }
1036
1037       snew(xrefp,atoms->nr);
1038       if (xref1 != NULL) {
1039         for(i=0; (i<nfit); i++) {
1040           copy_rvec(xref1[i],xrefp[ifit[i]]);
1041         }
1042       } else {
1043         /* The top coordinates are the fitting reference */
1044         for(i=0; (i<nfit); i++) {
1045           copy_rvec(xtop[ifit[i]],xrefp[ifit[i]]);
1046         }
1047         reset_x(nfit,ifit,atoms->nr,NULL,xrefp,w_rls);
1048       }
1049     }
1050     gmx_rmpbc_done(gpbc);
1051   }
1052
1053   if (bIndex) {
1054     printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n",natoms);
1055     get_index(atoms,indexfile,1,&i,&index,&grpname);
1056     if (i!=natoms)
1057       gmx_fatal(FARGS,"you selected a group with %d elements instead of %d",i,natoms);
1058       printf("\n");
1059   }
1060
1061   snew(sqrtm,natoms);
1062   if (bM && bDMA1) {
1063     proj_unit="u\\S1/2\\Nnm";
1064     for(i=0; (i<natoms); i++)
1065       sqrtm[i]=sqrt(atoms->atom[index[i]].m);
1066   }
1067   else {
1068     proj_unit="nm";
1069     for(i=0; (i<natoms); i++)
1070       sqrtm[i]=1.0;
1071   }
1072
1073   if (bVec2) {
1074     t=0;
1075     totmass=0;
1076     for(i=0; (i<natoms); i++)
1077       for(d=0;(d<DIM);d++) {
1078         t+=sqr((xav1[i][d]-xav2[i][d])*sqrtm[i]);
1079         totmass+=sqr(sqrtm[i]);
1080       }
1081     fprintf(stdout,"RMSD (without fit) between the two average structures:"
1082             " %.3f (nm)\n\n",sqrt(t/totmass));
1083   }
1084
1085   if (last==-1)
1086     last=natoms*DIM;
1087   if (first>-1) {
1088     if (bFirstToLast) {
1089       /* make an index from first to last */
1090       nout=last-first+1;
1091       snew(iout,nout);
1092       for(i=0; i<nout; i++)
1093         iout[i]=first-1+i;
1094     }
1095     else if (ThreeDPlotFile) {
1096       /* make an index of first+(0,1,2) and last */
1097       nout = bPDB3D ? 4 : 3;
1098       nout = min(last-first+1, nout);
1099       snew(iout,nout);
1100       iout[0]=first-1;
1101       iout[1]=first;
1102       if (nout>3)
1103         iout[2]=first+1;
1104       iout[nout-1]=last-1;
1105     }
1106     else {
1107       /* make an index of first and last */
1108       nout=2;
1109       snew(iout,nout);
1110       iout[0]=first-1;
1111       iout[1]=last-1;
1112     }
1113   }
1114   else {
1115     printf("Select eigenvectors for output, end your selection with 0\n");
1116     nout=-1;
1117     iout=NULL;
1118
1119     do {
1120       nout++;
1121       srenew(iout,nout+1);
1122       if(1 != scanf("%d",&iout[nout]))
1123       {
1124           gmx_fatal(FARGS,"Error reading user input");
1125       }
1126       iout[nout]--;
1127     }
1128     while (iout[nout]>=0);
1129
1130     printf("\n");
1131   }
1132   /* make an index of the eigenvectors which are present */
1133   snew(outvec,nout);
1134   noutvec=0;
1135   for(i=0; i<nout; i++)
1136     {
1137       j=0;
1138       while ((j<nvec1) && (eignr1[j]!=iout[i]))
1139         j++;
1140       if ((j<nvec1) && (eignr1[j]==iout[i]))
1141         {
1142           outvec[noutvec]=j;
1143           noutvec++;
1144         }
1145     }
1146   fprintf(stderr,"%d eigenvectors selected for output",noutvec);
1147   if (noutvec <= 100)
1148     {
1149       fprintf(stderr,":");
1150       for(j=0; j<noutvec; j++)
1151         fprintf(stderr," %d",eignr1[outvec[j]]+1);
1152     }
1153   fprintf(stderr,"\n");
1154
1155   if (CompFile)
1156     components(CompFile,natoms,eignr1,eigvec1,noutvec,outvec,oenv);
1157
1158   if (RmsfFile)
1159     rmsf(RmsfFile,natoms,sqrtm,eignr1,eigvec1,noutvec,outvec,eigval1,
1160          neig1,oenv);
1161
1162   if (bProj)
1163     project(bTraj ? opt2fn("-f",NFILE,fnm) : NULL,
1164             bTop ? &top : NULL,ePBC,topbox,
1165             ProjOnVecFile,TwoDPlotFile,ThreeDPlotFile,FilterFile,skip,
1166             ExtremeFile,bFirstLastSet,max,nextr,atoms,natoms,index,
1167             bFit1,xrefp,nfit,ifit,w_rls,
1168             sqrtm,xav1,eignr1,eigvec1,noutvec,outvec,bSplit,
1169             oenv);
1170
1171   if (OverlapFile)
1172     overlap(OverlapFile,natoms,
1173             eigvec1,nvec2,eignr2,eigvec2,noutvec,outvec,oenv);
1174
1175   if (InpMatFile)
1176     inprod_matrix(InpMatFile,natoms,
1177                   nvec1,eignr1,eigvec1,nvec2,eignr2,eigvec2,
1178                   bFirstLastSet,noutvec,outvec);
1179
1180   if (bCompare)
1181     compare(natoms,nvec1,eigvec1,nvec2,eigvec2,eigval1,neig1,eigval2,neig2);
1182
1183
1184   if (!CompFile && !bProj && !OverlapFile && !InpMatFile &&
1185           !bCompare && !bEntropy)
1186   {
1187     fprintf(stderr,"\nIf you want some output,"
1188             " set one (or two or ...) of the output file options\n");
1189   }
1190
1191
1192   view_all(oenv,NFILE, fnm);
1193
1194   thanx(stdout);
1195
1196   return 0;
1197 }
1198