src/tools/gmx_anaeig.c

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