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