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