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