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