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