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   atom_id *all_at;
 413   matrix  box;
 414   rvec    *xread,*x;
 415   real    t,inp,**inprod=NULL,min=0,max=0;
 416   char    str[STRLEN],str2[STRLEN],**ylabel,*c;
 417   real    fact;
 418   gmx_rmpbc_t  gpbc=NULL;
 419
 420   snew(x,natoms);
 421
 422   if (bExtrAll)
 423     noutvec_extr=noutvec;
 424   else
 425     noutvec_extr=1;
 426
 427
 428   if (trajfile) {
 429     snew(inprod,noutvec+1);
 430
 431     if (filterfile) {
 432       fprintf(stderr,"Writing a filtered trajectory to %s using eigenvectors\n",
 433               filterfile);
 434       for(i=0; i<noutvec; i++)
 435         fprintf(stderr,"%d ",outvec[i]+1);
 436       fprintf(stderr,"\n");
 437       out=open_trx(filterfile,"w");
 438     }
 439     snew_size=0;
 440     nfr=0;
 441     nframes=0;
 442     nat=read_first_x(oenv,&status,trajfile,&t,&xread,box);
 443     if (nat>atoms->nr)
 444       gmx_fatal(FARGS,"the number of atoms in your trajectory (%d) is larger than the number of atoms in your structure file (%d)",nat,atoms->nr);
 445     snew(all_at,nat);
 446
 447     if (top)
 448       gpbc = gmx_rmpbc_init(&top->idef,ePBC,nat,box);
 449     if (top)
 450     gmx_rmpbc_done(gpbc);
 451
 452     for(i=0; i<nat; i++)
 453       all_at[i]=i;
 454     do {
 455       if (nfr % skip == 0) {
 456         if (top)
 457           gmx_rmpbc(gpbc,nat,box,xread);
 458         if (nframes>=snew_size) {
 459           snew_size+=100;
 460           for(i=0; i<noutvec+1; i++)
 461             srenew(inprod[i],snew_size);
 462         }
 463         inprod[noutvec][nframes]=t;
 464         /* calculate x: a fitted struture of the selected atoms */
 465         if (bFit) {
 466           reset_x(nfit,ifit,nat,NULL,xread,w_rls);
 467           do_fit(nat,w_rls,xref,xread);
 468         }
 469         for (i=0; i<natoms; i++)
 470           copy_rvec(xread[index[i]],x[i]);
 471
 472         for(v=0; v<noutvec; v++) {
 473           vec=outvec[v];
 474           /* calculate (mass-weighted) projection */
 475           inp=0;
 476           for (i=0; i<natoms; i++) {
 477             inp+=(eigvec[vec][i][0]*(x[i][0]-xav[i][0])+
 478             eigvec[vec][i][1]*(x[i][1]-xav[i][1])+
 479             eigvec[vec][i][2]*(x[i][2]-xav[i][2]))*sqrtm[i];
 480           }
 481           inprod[v][nframes]=inp;
 482         }
 483         if (filterfile) {
 484           for(i=0; i<natoms; i++)
 485             for(d=0; d<DIM; d++) {
 486               /* misuse xread for output */
 487               xread[index[i]][d] = xav[i][d];
 488               for(v=0; v<noutvec; v++)
 489                 xread[index[i]][d] +=
 490                   inprod[v][nframes]*eigvec[outvec[v]][i][d]/sqrtm[i];
 491             }
 492           write_trx(out,natoms,index,atoms,0,t,box,xread,NULL,NULL);
 493         }
 494         nframes++;
 495       }
 496       nfr++;
 497     } while (read_next_x(oenv,status,&t,nat,xread,box));
 498     close_trx(status);
 499      sfree(x);
 500      if (filterfile)
 501        close_trx(out);
 502   }
 503   else
 504     snew(xread,atoms->nr);
 505
 506   if (top)
 507     gmx_rmpbc_done(gpbc);
 508
 509
 510   if (projfile) {
 511     snew(ylabel,noutvec);
 512     for(v=0; v<noutvec; v++) {
 513       sprintf(str,"vec %d",eignr[outvec[v]]+1);
 514       ylabel[v]=strdup(str);
 515     }
 516     sprintf(str,"projection on eigenvectors (%s)",proj_unit);
 517     write_xvgr_graphs(projfile, noutvec, 1, str, NULL, output_env_get_xvgr_tlabel(oenv),
 518                       (const char **)ylabel,
 519                       nframes, inprod[noutvec], inprod, NULL,
 520                       output_env_get_time_factor(oenv), FALSE, bSplit,oenv);
 521   }
 522
 523   if (twodplotfile) {
 524     sprintf(str,"projection on eigenvector %d (%s)",
 525             eignr[outvec[0]]+1,proj_unit);
 526     sprintf(str2,"projection on eigenvector %d (%s)",
 527             eignr[outvec[noutvec-1]]+1,proj_unit);
 528     xvgrout=xvgropen(twodplotfile,"2D projection of trajectory",str,str2,
 529                      oenv);
 530     for(i=0; i<nframes; i++) {
 531       if ( bSplit && i>0 && abs(inprod[noutvec][i])<1e-5 )
 532         fprintf(xvgrout,"&\n");
 533       fprintf(xvgrout,"%10.5f %10.5f\n",inprod[0][i],inprod[noutvec-1][i]);
 534     }
 535     ffclose(xvgrout);
 536   }
 537
 538   if (threedplotfile) {
 539     t_atoms atoms;
 540     rvec    *x;
 541     real    *b=NULL;
 542     matrix  box;
 543     char    *resnm,*atnm, pdbform[STRLEN];
 544     gmx_bool    bPDB, b4D;
 545     FILE    *out;
 546
 547     if (noutvec < 3)
 548       gmx_fatal(FARGS,"You have selected less than 3 eigenvectors");
 549
 550     /* initialize */
 551     bPDB = fn2ftp(threedplotfile)==efPDB;
 552     clear_mat(box);
 553     box[XX][XX] = box[YY][YY] = box[ZZ][ZZ] = 1;
 554
 555     b4D = bPDB && (noutvec >= 4);
 556     if (b4D) {
 557       fprintf(stderr, "You have selected four or more eigenvectors:\n"
 558               "fourth eigenvector will be plotted "
 559               "in bfactor field of pdb file\n");
 560       sprintf(str,"4D proj. of traj. on eigenv. %d, %d, %d and %d",
 561               eignr[outvec[0]]+1,eignr[outvec[1]]+1,
 562               eignr[outvec[2]]+1,eignr[outvec[3]]+1);
 563     } else {
 564     sprintf(str,"3D proj. of traj. on eigenv. %d, %d and %d",
 565             eignr[outvec[0]]+1,eignr[outvec[1]]+1,eignr[outvec[2]]+1);
 566     }
 567     init_t_atoms(&atoms,nframes,FALSE);
 568     snew(x,nframes);
 569     snew(b,nframes);
 570     atnm=strdup("C");
 571     resnm=strdup("PRJ");
 572
 573     if(nframes>10000)
 574       fact=10000.0/nframes;
 575     else
 576       fact=1.0;
 577
 578     for(i=0; i<nframes; i++) {
 579       atoms.atomname[i] = &atnm;
 580       atoms.atom[i].resind = i;
 581       atoms.resinfo[i].name = &resnm;
 582       atoms.resinfo[i].nr   = ceil(i*fact);
 583       atoms.resinfo[i].ic   = ' ';
 584       x[i][XX]=inprod[0][i];
 585       x[i][YY]=inprod[1][i];
 586       x[i][ZZ]=inprod[2][i];
 587       if (b4D)
 588         b[i]  =inprod[3][i];
 589     }
 590     if ( ( b4D || bSplit ) && bPDB ) {
 591       strcpy(pdbform,pdbformat);
 592       strcat(pdbform,"%8.4f%8.4f\n");
 593
 594       out=ffopen(threedplotfile,"w");
 595       fprintf(out,"HEADER    %s\n",str);
 596       if ( b4D )
 597         fprintf(out,"REMARK    %s\n","fourth dimension plotted as B-factor");
 598       j=0;
 599       for(i=0; i<atoms.nr; i++) {
 600         if ( j>0 && bSplit && abs(inprod[noutvec][i])<1e-5 ) {
 601           fprintf(out,"TER\n");
 602           j=0;
 603         }
 604         fprintf(out,pdbform,"ATOM",i+1,"C","PRJ",' ',j+1,
 605                 PR_VEC(10*x[i]), 1.0, 10*b[i]);
 606         if (j>0)
 607           fprintf(out,"CONECT%5d%5d\n", i, i+1);
 608         j++;
 609       }
 610       fprintf(out,"TER\n");
 611       ffclose(out);
 612     } else
 613       write_sto_conf(threedplotfile,str,&atoms,x,NULL,ePBC,box);
 614     free_t_atoms(&atoms,FALSE);
 615   }
 616
 617   if (extremefile) {
 618     if (extreme==0) {
 619       fprintf(stderr,"%11s %17s %17s\n","eigenvector","Minimum","Maximum");
 620       fprintf(stderr,
 621               "%11s %10s %10s %10s %10s\n","","value","time","value","time");
 622       snew(imin,noutvec_extr);
 623       snew(imax,noutvec_extr);
 624       for(v=0; v<noutvec_extr; v++) {
 625         for(i=0; i<nframes; i++) {
 626           if (inprod[v][i]<inprod[v][imin[v]])
 627             imin[v]=i;
 628           if (inprod[v][i]>inprod[v][imax[v]])
 629             imax[v]=i;
 630         }
 631         min=inprod[v][imin[v]];
 632         max=inprod[v][imax[v]];
 633         fprintf(stderr,"%7d     %10.6f %10.1f %10.6f %10.1f\n",
 634                 eignr[outvec[v]]+1,
 635                 min,inprod[noutvec][imin[v]],max,inprod[noutvec][imax[v]]);
 636       }
 637     }
 638     else {
 639       min=-extreme;
 640       max=+extreme;
 641     }
 642     /* build format string for filename: */
 643     strcpy(str,extremefile);/* copy filename */
 644     c=strrchr(str,'.'); /* find where extention begins */
 645     strcpy(str2,c); /* get extention */
 646     sprintf(c,"%%d%s",str2); /* append '%s' and extention to filename */
 647     for(v=0; v<noutvec_extr; v++) {
 648       /* make filename using format string */
 649       if (noutvec_extr==1)
 650         strcpy(str2,extremefile);
 651       else
 652         sprintf(str2,str,eignr[outvec[v]]+1);
 653       fprintf(stderr,"Writing %d frames along eigenvector %d to %s\n",
 654               nextr,outvec[v]+1,str2);
 655       out=open_trx(str2,"w");
 656       for(frame=0; frame<nextr; frame++) {
 657         if ((extreme==0) && (nextr<=3))
 658           for(i=0; i<natoms; i++) {
 659             atoms->resinfo[atoms->atom[index[i]].resind].chainid = 'A' + frame;
 660           }
 661         for(i=0; i<natoms; i++)
 662           for(d=0; d<DIM; d++)
 663             xread[index[i]][d] =
 664               (xav[i][d] + (min*(nextr-frame-1)+max*frame)/(nextr-1)
 665               *eigvec[outvec[v]][i][d]/sqrtm[i]);
 666         write_trx(out,natoms,index,atoms,0,frame,topbox,xread,NULL,NULL);
 667       }
 668       close_trx(out);
 669     }
 670   }
 671   fprintf(stderr,"\n");
 672 }
 673
 674 static void components(const char *outfile,int natoms,
 675                        int *eignr,rvec **eigvec,
 676                        int noutvec,int *outvec,
 677                        const output_env_t oenv)
 678 {
 679   int g,s,v,i;
 680   real *x,***y;
 681   char str[STRLEN],**ylabel;
 682
 683   fprintf(stderr,"Writing eigenvector components to %s\n",outfile);
 684
 685   snew(ylabel,noutvec);
 686   snew(y,noutvec);
 687   snew(x,natoms);
 688   for(i=0; i<natoms; i++)
 689     x[i]=i+1;
 690   for(g=0; g<noutvec; g++) {
 691     v=outvec[g];
 692     sprintf(str,"vec %d",eignr[v]+1);
 693     ylabel[g]=strdup(str);
 694     snew(y[g],4);
 695     for(s=0; s<4; s++) {
 696       snew(y[g][s],natoms);
 697     }
 698     for(i=0; i<natoms; i++) {
 699       y[g][0][i] = norm(eigvec[v][i]);
 700       for(s=0; s<3; s++)
 701         y[g][s+1][i] = eigvec[v][i][s];
 702     }
 703   }
 704   write_xvgr_graphs(outfile,noutvec,4,"Eigenvector components",
 705                     "black: total, red: x, green: y, blue: z",
 706                     "Atom number",(const char **)ylabel,
 707                     natoms,x,NULL,y,1,FALSE,FALSE,oenv);
 708   fprintf(stderr,"\n");
 709 }
 710
 711 static void rmsf(const char *outfile,int natoms,real *sqrtm,
 712                  int *eignr,rvec **eigvec,
 713                  int noutvec,int *outvec,
 714                  real *eigval, int neig,
 715                  const output_env_t oenv)
 716 {
 717   int   nrow,g,v,i;
 718   real  *x,**y;
 719   char  str[STRLEN],**ylabel;
 720
 721   for(i=0; i<neig; i++)
 722       if (eigval[i] < 0)
 723           eigval[i] = 0;
 724
 725   fprintf(stderr,"Writing rmsf to %s\n",outfile);
 726
 727   snew(ylabel,noutvec);
 728   snew(y,noutvec);
 729   snew(x,natoms);
 730   for(i=0; i<natoms; i++)
 731       x[i]=i+1;
 732   for(g=0; g<noutvec; g++) {
 733       v=outvec[g];
 734       if (eignr[v] >= neig)
 735           gmx_fatal(FARGS,"Selected vector %d is larger than the number of eigenvalues (%d)",eignr[v]+1,neig);
 736       sprintf(str,"vec %d",eignr[v]+1);
 737       ylabel[g]=strdup(str);
 738       snew(y[g],natoms);
 739       for(i=0; i<natoms; i++)
 740           y[g][i] = sqrt(eigval[eignr[v]]*norm2(eigvec[v][i]))/sqrtm[i];
 741   }
 742   write_xvgr_graphs(outfile,noutvec,1,"RMS fluctuation (nm) ",NULL,
 743                     "Atom number",(const char **)ylabel,
 744                     natoms,x,y,NULL,1,TRUE,FALSE,oenv);
 745   fprintf(stderr,"\n");
 746 }
 747
 748 int gmx_anaeig(int argc,char *argv[])
 749 {
 750   static const char *desc[] = {
 751     "[TT]g_anaeig[tt] analyzes eigenvectors. The eigenvectors can be of a",
 752     "covariance matrix ([TT]g_covar[tt]) or of a Normal Modes analysis",
 753     "([TT]g_nmeig[tt]).[PAR]",
 754
 755     "When a trajectory is projected on eigenvectors, all structures are",
 756     "fitted to the structure in the eigenvector file, if present, otherwise",
 757     "to the structure in the structure file. When no run input file is",
 758     "supplied, periodicity will not be taken into account. Most analyses",
 759     "are performed on eigenvectors [TT]-first[tt] to [TT]-last[tt], but when",
 760     "[TT]-first[tt] is set to -1 you will be prompted for a selection.[PAR]",
 761
 762     "[TT]-comp[tt]: plot the vector components per atom of eigenvectors",
 763     "[TT]-first[tt] to [TT]-last[tt].[PAR]",
 764
 765     "[TT]-rmsf[tt]: plot the RMS fluctuation per atom of eigenvectors",
 766     "[TT]-first[tt] to [TT]-last[tt] (requires [TT]-eig[tt]).[PAR]",
 767
 768     "[TT]-proj[tt]: calculate projections of a trajectory on eigenvectors",
 769     "[TT]-first[tt] to [TT]-last[tt].",
 770     "The projections of a trajectory on the eigenvectors of its",
 771     "covariance matrix are called principal components (pc's).",
 772     "It is often useful to check the cosine content of the pc's,",
 773     "since the pc's of random diffusion are cosines with the number",
 774     "of periods equal to half the pc index.",
 775     "The cosine content of the pc's can be calculated with the program",
 776     "[TT]g_analyze[tt].[PAR]",
 777
 778     "[TT]-2d[tt]: calculate a 2d projection of a trajectory on eigenvectors",
 779     "[TT]-first[tt] and [TT]-last[tt].[PAR]",
 780
 781     "[TT]-3d[tt]: calculate a 3d projection of a trajectory on the first",
 782     "three selected eigenvectors.[PAR]",
 783
 784     "[TT]-filt[tt]: filter the trajectory to show only the motion along",
 785     "eigenvectors [TT]-first[tt] to [TT]-last[tt].[PAR]",
 786
 787     "[TT]-extr[tt]: calculate the two extreme projections along a trajectory",
 788     "on the average structure and interpolate [TT]-nframes[tt] frames",
 789     "between them, or set your own extremes with [TT]-max[tt]. The",
 790     "eigenvector [TT]-first[tt] will be written unless [TT]-first[tt] and",
 791     "[TT]-last[tt] have been set explicitly, in which case all eigenvectors",
 792     "will be written to separate files. Chain identifiers will be added",
 793     "when writing a [TT].pdb[tt] file with two or three structures (you",
 794     "can use [TT]rasmol -nmrpdb[tt] to view such a [TT].pdb[tt] file).[PAR]",
 795
 796     "  Overlap calculations between covariance analysis:[BR]",
 797     "  NOTE: the analysis should use the same fitting structure[PAR]",
 798
 799     "[TT]-over[tt]: calculate the subspace overlap of the eigenvectors in",
 800     "file [TT]-v2[tt] with eigenvectors [TT]-first[tt] to [TT]-last[tt]",
 801     "in file [TT]-v[tt].[PAR]",
 802
 803     "[TT]-inpr[tt]: calculate a matrix of inner-products between",
 804     "eigenvectors in files [TT]-v[tt] and [TT]-v2[tt]. All eigenvectors",
 805     "of both files will be used unless [TT]-first[tt] and [TT]-last[tt]",
 806     "have been set explicitly.[PAR]",
 807
 808     "When [TT]-v[tt], [TT]-eig[tt], [TT]-v2[tt] and [TT]-eig2[tt] are given,",
 809     "a single number for the overlap between the covariance matrices is",
 810     "generated. The formulas are:[BR]",
 811     "        difference = sqrt(tr((sqrt(M1) - sqrt(M2))^2))[BR]",
 812     "normalized overlap = 1 - difference/sqrt(tr(M1) + tr(M2))[BR]",
 813     "     shape overlap = 1 - sqrt(tr((sqrt(M1/tr(M1)) - sqrt(M2/tr(M2)))^2))[BR]",
 814     "where M1 and M2 are the two covariance matrices and tr is the trace",
 815     "of a matrix. The numbers are proportional to the overlap of the square",
 816     "root of the fluctuations. The normalized overlap is the most useful",
 817     "number, it is 1 for identical matrices and 0 when the sampled",
 818     "subspaces are orthogonal.[PAR]",
 819     "When the [TT]-entropy[tt] flag is given an entropy estimate will be",
 820     "computed based on the Quasiharmonic approach and based on",
 821     "Schlitter's formula."
 822   };
 823   static int  first=1,last=8,skip=1,nextr=2,nskip=6;
 824   static real max=0.0,temp=298.15;
 825   static gmx_bool bSplit=FALSE,bEntropy=FALSE;
 826   t_pargs pa[] = {
 827     { "-first", FALSE, etINT, {&first},
 828       "First eigenvector for analysis (-1 is select)" },
 829     { "-last",  FALSE, etINT, {&last},
 830       "Last eigenvector for analysis (-1 is till the last)" },
 831     { "-skip",  FALSE, etINT, {&skip},
 832       "Only analyse every nr-th frame" },
 833     { "-max",  FALSE, etREAL, {&max},
 834       "Maximum for projection of the eigenvector on the average structure, "
 835       "max=0 gives the extremes" },
 836     { "-nframes",  FALSE, etINT, {&nextr},
 837       "Number of frames for the extremes output" },
 838     { "-split",   FALSE, etBOOL, {&bSplit},
 839       "Split eigenvector projections where time is zero" },
 840     { "-entropy", FALSE, etBOOL, {&bEntropy},
 841       "Compute entropy according to the Quasiharmonic formula or Schlitter's method." },
 842     { "-temp",    FALSE, etREAL, {&temp},
 843       "Temperature for entropy calculations" },
 844     { "-nevskip", FALSE, etINT, {&nskip},
 845       "Number of eigenvalues to skip when computing the entropy due to the quasi harmonic approximation. When you do a rotational and/or translational fit prior to the covariance analysis, you get 3 or 6 eigenvalues that are very close to zero, and which should not be taken into account when computing the entropy." }
 846   };
 847 #define NPA asize(pa)
 848
 849   FILE       *out;
 850   int        status,trjout;
 851   t_topology top;
 852   int        ePBC=-1;
 853   t_atoms    *atoms=NULL;
 854   rvec       *xtop,*xref1,*xref2,*xrefp=NULL;
 855   gmx_bool       bDMR1,bDMA1,bDMR2,bDMA2;
 856   int        nvec1,nvec2,*eignr1=NULL,*eignr2=NULL;
 857   rvec       *x,*xread,*xav1,*xav2,**eigvec1=NULL,**eigvec2=NULL;
 858   matrix     topbox;
 859   real       xid,totmass,*sqrtm,*w_rls,t,lambda;
 860   int        natoms,step;
 861   char       *grpname;
 862   const char *indexfile;
 863   char       title[STRLEN];
 864   int        i,j,d;
 865   int        nout,*iout,noutvec,*outvec,nfit;
 866   atom_id    *index,*ifit;
 867   const char *VecFile,*Vec2File,*topfile;
 868   const char *EigFile,*Eig2File;
 869   const char *CompFile,*RmsfFile,*ProjOnVecFile;
 870   const char *TwoDPlotFile,*ThreeDPlotFile;
 871   const char *FilterFile,*ExtremeFile;
 872   const char *OverlapFile,*InpMatFile;
 873   gmx_bool       bFit1,bFit2,bM,bIndex,bTPS,bTop,bVec2,bProj;
 874   gmx_bool       bFirstToLast,bFirstLastSet,bTraj,bCompare,bPDB3D;
 875   real       *eigval1=NULL,*eigval2=NULL;
 876   int        neig1,neig2;
 877   double     **xvgdata;
 878   output_env_t oenv;
 879   gmx_rmpbc_t  gpbc;
 880
 881   t_filenm fnm[] = {
 882     { efTRN, "-v",    "eigenvec",    ffREAD  },
 883     { efTRN, "-v2",   "eigenvec2",   ffOPTRD },
 884     { efTRX, "-f",    NULL,          ffOPTRD },
 885     { efTPS, NULL,    NULL,          ffOPTRD },
 886     { efNDX, NULL,    NULL,          ffOPTRD },
 887     { efXVG, "-eig", "eigenval",     ffOPTRD },
 888     { efXVG, "-eig2", "eigenval2",   ffOPTRD },
 889     { efXVG, "-comp", "eigcomp",     ffOPTWR },
 890     { efXVG, "-rmsf", "eigrmsf",     ffOPTWR },
 891     { efXVG, "-proj", "proj",        ffOPTWR },
 892     { efXVG, "-2d",   "2dproj",      ffOPTWR },
 893     { efSTO, "-3d",   "3dproj.pdb",  ffOPTWR },
 894     { efTRX, "-filt", "filtered",    ffOPTWR },
 895     { efTRX, "-extr", "extreme.pdb", ffOPTWR },
 896     { efXVG, "-over", "overlap",     ffOPTWR },
 897     { efXPM, "-inpr", "inprod",      ffOPTWR }
 898   };
 899 #define NFILE asize(fnm)
 900
 901   CopyRight(stderr,argv[0]);
 902   parse_common_args(&argc,argv,
 903                     PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW | PCA_BE_NICE ,
 904                     NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
 905
 906   indexfile=ftp2fn_null(efNDX,NFILE,fnm);
 907
 908   VecFile         = opt2fn("-v",NFILE,fnm);
 909   Vec2File        = opt2fn_null("-v2",NFILE,fnm);
 910   topfile         = ftp2fn(efTPS,NFILE,fnm);
 911   EigFile         = opt2fn_null("-eig",NFILE,fnm);
 912   Eig2File        = opt2fn_null("-eig2",NFILE,fnm);
 913   CompFile        = opt2fn_null("-comp",NFILE,fnm);
 914   RmsfFile        = opt2fn_null("-rmsf",NFILE,fnm);
 915   ProjOnVecFile   = opt2fn_null("-proj",NFILE,fnm);
 916   TwoDPlotFile    = opt2fn_null("-2d",NFILE,fnm);
 917   ThreeDPlotFile  = opt2fn_null("-3d",NFILE,fnm);
 918   FilterFile      = opt2fn_null("-filt",NFILE,fnm);
 919   ExtremeFile     = opt2fn_null("-extr",NFILE,fnm);
 920   OverlapFile     = opt2fn_null("-over",NFILE,fnm);
 921   InpMatFile      = ftp2fn_null(efXPM,NFILE,fnm);
 922
 923   bTop   = fn2bTPX(topfile);
 924   bProj  = ProjOnVecFile || TwoDPlotFile || ThreeDPlotFile
 925     || FilterFile || ExtremeFile;
 926   bFirstLastSet  =
 927     opt2parg_bSet("-first",NPA,pa) && opt2parg_bSet("-last",NPA,pa);
 928   bFirstToLast = CompFile || RmsfFile || ProjOnVecFile || FilterFile ||
 929     OverlapFile || ((ExtremeFile || InpMatFile) && bFirstLastSet);
 930   bVec2  = Vec2File || OverlapFile || InpMatFile;
 931   bM     = RmsfFile || bProj;
 932   bTraj  = ProjOnVecFile || FilterFile || (ExtremeFile && (max==0))
 933     || TwoDPlotFile || ThreeDPlotFile;
 934   bIndex = bM || bProj;
 935   bTPS   = ftp2bSet(efTPS,NFILE,fnm) || bM || bTraj ||
 936     FilterFile  || (bIndex && indexfile);
 937   bCompare = Vec2File || Eig2File;
 938   bPDB3D = fn2ftp(ThreeDPlotFile)==efPDB;
 939
 940   read_eigenvectors(VecFile,&natoms,&bFit1,
 941                     &xref1,&bDMR1,&xav1,&bDMA1,
 942                     &nvec1,&eignr1,&eigvec1,&eigval1);
 943   neig1 = DIM*natoms;
 944
 945   /* Overwrite eigenvalues from separate files if the user provides them */
 946   if (EigFile != NULL) {
 947     int neig_tmp = read_xvg(EigFile,&xvgdata,&i);
 948     if (neig_tmp != neig1)
 949       fprintf(stderr,"Warning: number of eigenvalues in xvg file (%d) does not mtch trr file (%d)\n",neig1,natoms);
 950     neig1 = neig_tmp;
 951     srenew(eigval1,neig1);
 952     for(j=0;j<neig1;j++) {
 953       real tmp = eigval1[j];
 954       eigval1[j]=xvgdata[1][j];
 955       if (debug && (eigval1[j] != tmp))
 956         fprintf(debug,"Replacing eigenvalue %d. From trr: %10g, from xvg: %10g\n",
 957                 j,tmp,eigval1[j]);
 958     }
 959     for(j=0;j<i;j++)
 960       sfree(xvgdata[j]);
 961     sfree(xvgdata);
 962     fprintf(stderr,"Read %d eigenvalues from %s\n",neig1,EigFile);
 963   }
 964
 965   if (bEntropy) {
 966     calc_entropy_qh(stdout,neig1,eigval1,temp,nskip);
 967     calc_entropy_schlitter(stdout,neig1,nskip,eigval1,temp);
 968   }
 969
 970   if (bVec2) {
 971     if (!Vec2File)
 972       gmx_fatal(FARGS,"Need a second eigenvector file to do this analysis.");
 973     read_eigenvectors(Vec2File,&neig2,&bFit2,
 974                       &xref2,&bDMR2,&xav2,&bDMA2,&nvec2,&eignr2,&eigvec2,&eigval2);
 975
 976     neig2 = DIM*neig2;
 977     if (neig2 != neig1)
 978       gmx_fatal(FARGS,"Dimensions in the eigenvector files don't match");
 979   }
 980
 981   if(Eig2File != NULL) {
 982     neig2 = read_xvg(Eig2File,&xvgdata,&i);
 983     srenew(eigval2,neig2);
 984     for(j=0;j<neig2;j++)
 985       eigval2[j]=xvgdata[1][j];
 986     for(j=0;j<i;j++)
 987       sfree(xvgdata[j]);
 988     sfree(xvgdata);
 989     fprintf(stderr,"Read %d eigenvalues from %s\n",neig2,Eig2File);
 990   }
 991
 992
 993   if ((!bFit1 || xref1) && !bDMR1 && !bDMA1)
 994     bM=FALSE;
 995   if ((xref1==NULL) && (bM || bTraj))
 996     bTPS=TRUE;
 997
 998   xtop=NULL;
 999   nfit=0;
1000   ifit=NULL;
1001   w_rls=NULL;
1002
1003   if (!bTPS) {
1004     bTop=FALSE;
1005   } else {
1006     bTop=read_tps_conf(ftp2fn(efTPS,NFILE,fnm),
1007                        title,&top,&ePBC,&xtop,NULL,topbox,bM);
1008     atoms=&top.atoms;
1009     gpbc = gmx_rmpbc_init(&top.idef,ePBC,atoms->nr,topbox);
1010     gmx_rmpbc(gpbc,atoms->nr,topbox,xtop);
1011     /* Fitting is only required for the projection */
1012     if (bProj && bFit1) {
1013       if (xref1 == NULL) {
1014           printf("\nNote: the structure in %s should be the same\n"
1015                  "      as the one used for the fit in g_covar\n",topfile);
1016       }
1017       printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
1018       get_index(atoms,indexfile,1,&nfit,&ifit,&grpname);
1019
1020       snew(w_rls,atoms->nr);
1021       for(i=0; (i<nfit); i++) {
1022         if (bDMR1) {
1023           w_rls[ifit[i]] = atoms->atom[ifit[i]].m;
1024         } else {
1025           w_rls[ifit[i]] = 1.0;
1026         }
1027       }
1028
1029       snew(xrefp,atoms->nr);
1030       if (xref1 != NULL) {
1031         for(i=0; (i<nfit); i++) {
1032           copy_rvec(xref1[i],xrefp[ifit[i]]);
1033         }
1034       } else {
1035         /* The top coordinates are the fitting reference */
1036         for(i=0; (i<nfit); i++) {
1037           copy_rvec(xtop[ifit[i]],xrefp[ifit[i]]);
1038         }
1039         reset_x(nfit,ifit,atoms->nr,NULL,xrefp,w_rls);
1040       }
1041     }
1042     gmx_rmpbc_done(gpbc);
1043   }
1044
1045   if (bIndex) {
1046     printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n",natoms);
1047     get_index(atoms,indexfile,1,&i,&index,&grpname);
1048     if (i!=natoms)
1049       gmx_fatal(FARGS,"you selected a group with %d elements instead of %d",i,natoms);
1050       printf("\n");
1051   }
1052
1053   snew(sqrtm,natoms);
1054   if (bM && bDMA1) {
1055     proj_unit="u\\S1/2\\Nnm";
1056     for(i=0; (i<natoms); i++)
1057       sqrtm[i]=sqrt(atoms->atom[index[i]].m);
1058   }
1059   else {
1060     proj_unit="nm";
1061     for(i=0; (i<natoms); i++)
1062       sqrtm[i]=1.0;
1063   }
1064
1065   if (bVec2) {
1066     t=0;
1067     totmass=0;
1068     for(i=0; (i<natoms); i++)
1069       for(d=0;(d<DIM);d++) {
1070         t+=sqr((xav1[i][d]-xav2[i][d])*sqrtm[i]);
1071         totmass+=sqr(sqrtm[i]);
1072       }
1073     fprintf(stdout,"RMSD (without fit) between the two average structures:"
1074             " %.3f (nm)\n\n",sqrt(t/totmass));
1075   }
1076
1077   if (last==-1)
1078     last=natoms*DIM;
1079   if (first>-1) {
1080     if (bFirstToLast) {
1081       /* make an index from first to last */
1082       nout=last-first+1;
1083       snew(iout,nout);
1084       for(i=0; i<nout; i++)
1085         iout[i]=first-1+i;
1086     }
1087     else if (ThreeDPlotFile) {
1088       /* make an index of first+(0,1,2) and last */
1089       nout = bPDB3D ? 4 : 3;
1090       nout = min(last-first+1, nout);
1091       snew(iout,nout);
1092       iout[0]=first-1;
1093       iout[1]=first;
1094       if (nout>3)
1095         iout[2]=first+1;
1096       iout[nout-1]=last-1;
1097     }
1098     else {
1099       /* make an index of first and last */
1100       nout=2;
1101       snew(iout,nout);
1102       iout[0]=first-1;
1103       iout[1]=last-1;
1104     }
1105   }
1106   else {
1107     printf("Select eigenvectors for output, end your selection with 0\n");
1108     nout=-1;
1109     iout=NULL;
1110
1111     do {
1112       nout++;
1113       srenew(iout,nout+1);
1114       if(1 != scanf("%d",&iout[nout]))
1115       {
1116           gmx_fatal(FARGS,"Error reading user input");
1117       }
1118       iout[nout]--;
1119     }
1120     while (iout[nout]>=0);
1121
1122     printf("\n");
1123   }
1124   /* make an index of the eigenvectors which are present */
1125   snew(outvec,nout);
1126   noutvec=0;
1127   for(i=0; i<nout; i++)
1128     {
1129       j=0;
1130       while ((j<nvec1) && (eignr1[j]!=iout[i]))
1131         j++;
1132       if ((j<nvec1) && (eignr1[j]==iout[i]))
1133         {
1134           outvec[noutvec]=j;
1135           noutvec++;
1136         }
1137     }
1138   fprintf(stderr,"%d eigenvectors selected for output",noutvec);
1139   if (noutvec <= 100)
1140     {
1141       fprintf(stderr,":");
1142       for(j=0; j<noutvec; j++)
1143         fprintf(stderr," %d",eignr1[outvec[j]]+1);
1144     }
1145   fprintf(stderr,"\n");
1146
1147   if (CompFile)
1148     components(CompFile,natoms,eignr1,eigvec1,noutvec,outvec,oenv);
1149
1150   if (RmsfFile)
1151     rmsf(RmsfFile,natoms,sqrtm,eignr1,eigvec1,noutvec,outvec,eigval1,
1152          neig1,oenv);
1153
1154   if (bProj)
1155     project(bTraj ? opt2fn("-f",NFILE,fnm) : NULL,
1156             bTop ? &top : NULL,ePBC,topbox,
1157             ProjOnVecFile,TwoDPlotFile,ThreeDPlotFile,FilterFile,skip,
1158             ExtremeFile,bFirstLastSet,max,nextr,atoms,natoms,index,
1159             bFit1,xrefp,nfit,ifit,w_rls,
1160             sqrtm,xav1,eignr1,eigvec1,noutvec,outvec,bSplit,
1161             oenv);
1162
1163   if (OverlapFile)
1164     overlap(OverlapFile,natoms,
1165             eigvec1,nvec2,eignr2,eigvec2,noutvec,outvec,oenv);
1166
1167   if (InpMatFile)
1168     inprod_matrix(InpMatFile,natoms,
1169                   nvec1,eignr1,eigvec1,nvec2,eignr2,eigvec2,
1170                   bFirstLastSet,noutvec,outvec);
1171
1172   if (bCompare)
1173     compare(natoms,nvec1,eigvec1,nvec2,eigvec2,eigval1,neig1,eigval2,neig2);
1174
1175
1176   if (!CompFile && !bProj && !OverlapFile && !InpMatFile &&
1177           !bCompare && !bEntropy)
1178   {
1179     fprintf(stderr,"\nIf you want some output,"
1180             " set one (or two or ...) of the output file options\n");
1181   }
1182
1183
1184   view_all(oenv,NFILE, fnm);
1185
1186   thanx(stdout);
1187
1188   return 0;
1189 }
1190