src/gromacs/gmxana/gmx_covar.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 <string.h>
  43 #include <time.h>
  44
  45 #ifdef HAVE_SYS_TIME_H
  46 #include <sys/time.h>
  47 #endif
  48
  49 #include "gromacs/commandline/pargs.h"
  50 #include "sysstuff.h"
  51 #include "typedefs.h"
  52 #include "gromacs/utility/smalloc.h"
  53 #include "macros.h"
  54 #include "vec.h"
  55 #include "pbc.h"
  56 #include "gromacs/fileio/futil.h"
  57 #include "index.h"
  58 #include "gromacs/fileio/confio.h"
  59 #include "gromacs/fileio/trnio.h"
  60 #include "mshift.h"
  61 #include "xvgr.h"
  62 #include "rmpbc.h"
  63 #include "txtdump.h"
  64 #include "gromacs/fileio/matio.h"
  65 #include "eigio.h"
  66 #include "physics.h"
  67 #include "gmx_ana.h"
  68 #include "gromacs/utility/cstringutil.h"
  69 #include "gromacs/fileio/trxio.h"
  70
  71 #include "gromacs/linearalgebra/eigensolver.h"
  72 #include "gromacs/math/do_fit.h"
  73 #include "gromacs/legacyheaders/gmx_fatal.h"
  74
  75 int gmx_covar(int argc, char *argv[])
  76 {
  77     const char     *desc[] = {
  78         "[THISMODULE] calculates and diagonalizes the (mass-weighted)",
  79         "covariance matrix.",
  80         "All structures are fitted to the structure in the structure file.",
  81         "When this is not a run input file periodicity will not be taken into",
  82         "account. When the fit and analysis groups are identical and the analysis",
  83         "is non mass-weighted, the fit will also be non mass-weighted.",
  84         "[PAR]",
  85         "The eigenvectors are written to a trajectory file ([TT]-v[tt]).",
  86         "When the same atoms are used for the fit and the covariance analysis,",
  87         "the reference structure for the fit is written first with t=-1.",
  88         "The average (or reference when [TT]-ref[tt] is used) structure is",
  89         "written with t=0, the eigenvectors",
  90         "are written as frames with the eigenvector number as timestamp.",
  91         "[PAR]",
  92         "The eigenvectors can be analyzed with [gmx-anaeig].",
  93         "[PAR]",
  94         "Option [TT]-ascii[tt] writes the whole covariance matrix to",
  95         "an ASCII file. The order of the elements is: x1x1, x1y1, x1z1, x1x2, ...",
  96         "[PAR]",
  97         "Option [TT]-xpm[tt] writes the whole covariance matrix to an [TT].xpm[tt] file.",
  98         "[PAR]",
  99         "Option [TT]-xpma[tt] writes the atomic covariance matrix to an [TT].xpm[tt] file,",
 100         "i.e. for each atom pair the sum of the xx, yy and zz covariances is",
 101         "written.",
 102         "[PAR]",
 103         "Note that the diagonalization of a matrix requires memory and time",
 104         "that will increase at least as fast as than the square of the number",
 105         "of atoms involved. It is easy to run out of memory, in which",
 106         "case this tool will probably exit with a 'Segmentation fault'. You",
 107         "should consider carefully whether a reduced set of atoms will meet",
 108         "your needs for lower costs."
 109     };
 110     static gmx_bool bFit = TRUE, bRef = FALSE, bM = FALSE, bPBC = TRUE;
 111     static int      end  = -1;
 112     t_pargs         pa[] = {
 113         { "-fit",  FALSE, etBOOL, {&bFit},
 114           "Fit to a reference structure"},
 115         { "-ref",  FALSE, etBOOL, {&bRef},
 116           "Use the deviation from the conformation in the structure file instead of from the average" },
 117         { "-mwa",  FALSE, etBOOL, {&bM},
 118           "Mass-weighted covariance analysis"},
 119         { "-last",  FALSE, etINT, {&end},
 120           "Last eigenvector to write away (-1 is till the last)" },
 121         { "-pbc",  FALSE,  etBOOL, {&bPBC},
 122           "Apply corrections for periodic boundary conditions" }
 123     };
 124     FILE           *out = NULL; /* initialization makes all compilers happy */
 125     t_trxstatus    *status;
 126     t_trxstatus    *trjout;
 127     t_topology      top;
 128     int             ePBC;
 129     t_atoms        *atoms;
 130     rvec           *x, *xread, *xref, *xav, *xproj;
 131     matrix          box, zerobox;
 132     real           *sqrtm, *mat, *eigenvalues, sum, trace, inv_nframes;
 133     real            t, tstart, tend, **mat2;
 134     real            xj, *w_rls = NULL;
 135     real            min, max, *axis;
 136     int             ntopatoms, step;
 137     int             natoms, nat, count, nframes0, nframes, nlevels;
 138     gmx_int64_t     ndim, i, j, k, l;
 139     int             WriteXref;
 140     const char     *fitfile, *trxfile, *ndxfile;
 141     const char     *eigvalfile, *eigvecfile, *averfile, *logfile;
 142     const char     *asciifile, *xpmfile, *xpmafile;
 143     char            str[STRLEN], *fitname, *ananame, *pcwd;
 144     int             d, dj, nfit;
 145     atom_id        *index, *ifit;
 146     gmx_bool        bDiffMass1, bDiffMass2;
 147     time_t          now;
 148     char            timebuf[STRLEN];
 149     t_rgb           rlo, rmi, rhi;
 150     real           *eigenvectors;
 151     output_env_t    oenv;
 152     gmx_rmpbc_t     gpbc = NULL;
 153
 154     t_filenm        fnm[] = {
 155         { efTRX, "-f",  NULL, ffREAD },
 156         { efTPS, NULL,  NULL, ffREAD },
 157         { efNDX, NULL,  NULL, ffOPTRD },
 158         { efXVG, NULL,  "eigenval", ffWRITE },
 159         { efTRN, "-v",  "eigenvec", ffWRITE },
 160         { efSTO, "-av", "average.pdb", ffWRITE },
 161         { efLOG, NULL,  "covar", ffWRITE },
 162         { efDAT, "-ascii", "covar", ffOPTWR },
 163         { efXPM, "-xpm", "covar", ffOPTWR },
 164         { efXPM, "-xpma", "covara", ffOPTWR }
 165     };
 166 #define NFILE asize(fnm)
 167
 168     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
 169                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
 170     {
 171         return 0;
 172     }
 173
 174     clear_mat(zerobox);
 175
 176     fitfile    = ftp2fn(efTPS, NFILE, fnm);
 177     trxfile    = ftp2fn(efTRX, NFILE, fnm);
 178     ndxfile    = ftp2fn_null(efNDX, NFILE, fnm);
 179     eigvalfile = ftp2fn(efXVG, NFILE, fnm);
 180     eigvecfile = ftp2fn(efTRN, NFILE, fnm);
 181     averfile   = ftp2fn(efSTO, NFILE, fnm);
 182     logfile    = ftp2fn(efLOG, NFILE, fnm);
 183     asciifile  = opt2fn_null("-ascii", NFILE, fnm);
 184     xpmfile    = opt2fn_null("-xpm", NFILE, fnm);
 185     xpmafile   = opt2fn_null("-xpma", NFILE, fnm);
 186
 187     read_tps_conf(fitfile, str, &top, &ePBC, &xref, NULL, box, TRUE);
 188     atoms = &top.atoms;
 189
 190     if (bFit)
 191     {
 192         printf("\nChoose a group for the least squares fit\n");
 193         get_index(atoms, ndxfile, 1, &nfit, &ifit, &fitname);
 194         if (nfit < 3)
 195         {
 196             gmx_fatal(FARGS, "Need >= 3 points to fit!\n");
 197         }
 198     }
 199     else
 200     {
 201         nfit = 0;
 202     }
 203     printf("\nChoose a group for the covariance analysis\n");
 204     get_index(atoms, ndxfile, 1, &natoms, &index, &ananame);
 205
 206     bDiffMass1 = FALSE;
 207     if (bFit)
 208     {
 209         snew(w_rls, atoms->nr);
 210         for (i = 0; (i < nfit); i++)
 211         {
 212             w_rls[ifit[i]] = atoms->atom[ifit[i]].m;
 213             if (i)
 214             {
 215                 bDiffMass1 = bDiffMass1 || (w_rls[ifit[i]] != w_rls[ifit[i-1]]);
 216             }
 217         }
 218     }
 219     bDiffMass2 = FALSE;
 220     snew(sqrtm, natoms);
 221     for (i = 0; (i < natoms); i++)
 222     {
 223         if (bM)
 224         {
 225             sqrtm[i] = sqrt(atoms->atom[index[i]].m);
 226             if (i)
 227             {
 228                 bDiffMass2 = bDiffMass2 || (sqrtm[i] != sqrtm[i-1]);
 229             }
 230         }
 231         else
 232         {
 233             sqrtm[i] = 1.0;
 234         }
 235     }
 236
 237     if (bFit && bDiffMass1 && !bDiffMass2)
 238     {
 239         bDiffMass1 = natoms != nfit;
 240         i          = 0;
 241         for (i = 0; (i < natoms) && !bDiffMass1; i++)
 242         {
 243             bDiffMass1 = index[i] != ifit[i];
 244         }
 245         if (!bDiffMass1)
 246         {
 247             fprintf(stderr, "\n"
 248                     "Note: the fit and analysis group are identical,\n"
 249                     "      while the fit is mass weighted and the analysis is not.\n"
 250                     "      Making the fit non mass weighted.\n\n");
 251             for (i = 0; (i < nfit); i++)
 252             {
 253                 w_rls[ifit[i]] = 1.0;
 254             }
 255         }
 256     }
 257
 258     /* Prepare reference frame */
 259     if (bPBC)
 260     {
 261         gpbc = gmx_rmpbc_init(&top.idef, ePBC, atoms->nr);
 262         gmx_rmpbc(gpbc, atoms->nr, box, xref);
 263     }
 264     if (bFit)
 265     {
 266         reset_x(nfit, ifit, atoms->nr, NULL, xref, w_rls);
 267     }
 268
 269     snew(x, natoms);
 270     snew(xav, natoms);
 271     ndim = natoms*DIM;
 272     if (sqrt(GMX_INT64_MAX) < ndim)
 273     {
 274         gmx_fatal(FARGS, "Number of degrees of freedoms to large for matrix.\n");
 275     }
 276     snew(mat, ndim*ndim);
 277
 278     fprintf(stderr, "Calculating the average structure ...\n");
 279     nframes0 = 0;
 280     nat      = read_first_x(oenv, &status, trxfile, &t, &xread, box);
 281     if (nat != atoms->nr)
 282     {
 283         fprintf(stderr, "\nWARNING: number of atoms in tpx (%d) and trajectory (%d) do not match\n", natoms, nat);
 284     }
 285     do
 286     {
 287         nframes0++;
 288         /* calculate x: a fitted struture of the selected atoms */
 289         if (bPBC)
 290         {
 291             gmx_rmpbc(gpbc, nat, box, xread);
 292         }
 293         if (bFit)
 294         {
 295             reset_x(nfit, ifit, nat, NULL, xread, w_rls);
 296             do_fit(nat, w_rls, xref, xread);
 297         }
 298         for (i = 0; i < natoms; i++)
 299         {
 300             rvec_inc(xav[i], xread[index[i]]);
 301         }
 302     }
 303     while (read_next_x(oenv, status, &t, xread, box));
 304     close_trj(status);
 305
 306     inv_nframes = 1.0/nframes0;
 307     for (i = 0; i < natoms; i++)
 308     {
 309         for (d = 0; d < DIM; d++)
 310         {
 311             xav[i][d]         *= inv_nframes;
 312             xread[index[i]][d] = xav[i][d];
 313         }
 314     }
 315     write_sto_conf_indexed(opt2fn("-av", NFILE, fnm), "Average structure",
 316                            atoms, xread, NULL, epbcNONE, zerobox, natoms, index);
 317     sfree(xread);
 318
 319     fprintf(stderr, "Constructing covariance matrix (%dx%d) ...\n", (int)ndim, (int)ndim);
 320     nframes = 0;
 321     nat     = read_first_x(oenv, &status, trxfile, &t, &xread, box);
 322     tstart  = t;
 323     do
 324     {
 325         nframes++;
 326         tend = t;
 327         /* calculate x: a (fitted) structure of the selected atoms */
 328         if (bPBC)
 329         {
 330             gmx_rmpbc(gpbc, nat, box, xread);
 331         }
 332         if (bFit)
 333         {
 334             reset_x(nfit, ifit, nat, NULL, xread, w_rls);
 335             do_fit(nat, w_rls, xref, xread);
 336         }
 337         if (bRef)
 338         {
 339             for (i = 0; i < natoms; i++)
 340             {
 341                 rvec_sub(xread[index[i]], xref[index[i]], x[i]);
 342             }
 343         }
 344         else
 345         {
 346             for (i = 0; i < natoms; i++)
 347             {
 348                 rvec_sub(xread[index[i]], xav[i], x[i]);
 349             }
 350         }
 351
 352         for (j = 0; j < natoms; j++)
 353         {
 354             for (dj = 0; dj < DIM; dj++)
 355             {
 356                 k  = ndim*(DIM*j+dj);
 357                 xj = x[j][dj];
 358                 for (i = j; i < natoms; i++)
 359                 {
 360                     l = k+DIM*i;
 361                     for (d = 0; d < DIM; d++)
 362                     {
 363                         mat[l+d] += x[i][d]*xj;
 364                     }
 365                 }
 366             }
 367         }
 368     }
 369     while (read_next_x(oenv, status, &t, xread, box) &&
 370            (bRef || nframes < nframes0));
 371     close_trj(status);
 372     gmx_rmpbc_done(gpbc);
 373
 374     fprintf(stderr, "Read %d frames\n", nframes);
 375
 376     if (bRef)
 377     {
 378         /* copy the reference structure to the ouput array x */
 379         snew(xproj, natoms);
 380         for (i = 0; i < natoms; i++)
 381         {
 382             copy_rvec(xref[index[i]], xproj[i]);
 383         }
 384     }
 385     else
 386     {
 387         xproj = xav;
 388     }
 389
 390     /* correct the covariance matrix for the mass */
 391     inv_nframes = 1.0/nframes;
 392     for (j = 0; j < natoms; j++)
 393     {
 394         for (dj = 0; dj < DIM; dj++)
 395         {
 396             for (i = j; i < natoms; i++)
 397             {
 398                 k = ndim*(DIM*j+dj)+DIM*i;
 399                 for (d = 0; d < DIM; d++)
 400                 {
 401                     mat[k+d] = mat[k+d]*inv_nframes*sqrtm[i]*sqrtm[j];
 402                 }
 403             }
 404         }
 405     }
 406
 407     /* symmetrize the matrix */
 408     for (j = 0; j < ndim; j++)
 409     {
 410         for (i = j; i < ndim; i++)
 411         {
 412             mat[ndim*i+j] = mat[ndim*j+i];
 413         }
 414     }
 415
 416     trace = 0;
 417     for (i = 0; i < ndim; i++)
 418     {
 419         trace += mat[i*ndim+i];
 420     }
 421     fprintf(stderr, "\nTrace of the covariance matrix: %g (%snm^2)\n",
 422             trace, bM ? "u " : "");
 423
 424     if (asciifile)
 425     {
 426         out = gmx_ffopen(asciifile, "w");
 427         for (j = 0; j < ndim; j++)
 428         {
 429             for (i = 0; i < ndim; i += 3)
 430             {
 431                 fprintf(out, "%g %g %g\n",
 432                         mat[ndim*j+i], mat[ndim*j+i+1], mat[ndim*j+i+2]);
 433             }
 434         }
 435         gmx_ffclose(out);
 436     }
 437
 438     if (xpmfile)
 439     {
 440         min = 0;
 441         max = 0;
 442         snew(mat2, ndim);
 443         for (j = 0; j < ndim; j++)
 444         {
 445             mat2[j] = &(mat[ndim*j]);
 446             for (i = 0; i <= j; i++)
 447             {
 448                 if (mat2[j][i] < min)
 449                 {
 450                     min = mat2[j][i];
 451                 }
 452                 if (mat2[j][j] > max)
 453                 {
 454                     max = mat2[j][i];
 455                 }
 456             }
 457         }
 458         snew(axis, ndim);
 459         for (i = 0; i < ndim; i++)
 460         {
 461             axis[i] = i+1;
 462         }
 463         rlo.r   = 0; rlo.g = 0; rlo.b = 1;
 464         rmi.r   = 1; rmi.g = 1; rmi.b = 1;
 465         rhi.r   = 1; rhi.g = 0; rhi.b = 0;
 466         out     = gmx_ffopen(xpmfile, "w");
 467         nlevels = 80;
 468         write_xpm3(out, 0, "Covariance", bM ? "u nm^2" : "nm^2",
 469                    "dim", "dim", ndim, ndim, axis, axis,
 470                    mat2, min, 0.0, max, rlo, rmi, rhi, &nlevels);
 471         gmx_ffclose(out);
 472         sfree(axis);
 473         sfree(mat2);
 474     }
 475
 476     if (xpmafile)
 477     {
 478         min = 0;
 479         max = 0;
 480         snew(mat2, ndim/DIM);
 481         for (i = 0; i < ndim/DIM; i++)
 482         {
 483             snew(mat2[i], ndim/DIM);
 484         }
 485         for (j = 0; j < ndim/DIM; j++)
 486         {
 487             for (i = 0; i <= j; i++)
 488             {
 489                 mat2[j][i] = 0;
 490                 for (d = 0; d < DIM; d++)
 491                 {
 492                     mat2[j][i] += mat[ndim*(DIM*j+d)+DIM*i+d];
 493                 }
 494                 if (mat2[j][i] < min)
 495                 {
 496                     min = mat2[j][i];
 497                 }
 498                 if (mat2[j][j] > max)
 499                 {
 500                     max = mat2[j][i];
 501                 }
 502                 mat2[i][j] = mat2[j][i];
 503             }
 504         }
 505         snew(axis, ndim/DIM);
 506         for (i = 0; i < ndim/DIM; i++)
 507         {
 508             axis[i] = i+1;
 509         }
 510         rlo.r   = 0; rlo.g = 0; rlo.b = 1;
 511         rmi.r   = 1; rmi.g = 1; rmi.b = 1;
 512         rhi.r   = 1; rhi.g = 0; rhi.b = 0;
 513         out     = gmx_ffopen(xpmafile, "w");
 514         nlevels = 80;
 515         write_xpm3(out, 0, "Covariance", bM ? "u nm^2" : "nm^2",
 516                    "atom", "atom", ndim/DIM, ndim/DIM, axis, axis,
 517                    mat2, min, 0.0, max, rlo, rmi, rhi, &nlevels);
 518         gmx_ffclose(out);
 519         sfree(axis);
 520         for (i = 0; i < ndim/DIM; i++)
 521         {
 522             sfree(mat2[i]);
 523         }
 524         sfree(mat2);
 525     }
 526
 527
 528     /* call diagonalization routine */
 529
 530     snew(eigenvalues, ndim);
 531     snew(eigenvectors, ndim*ndim);
 532
 533     memcpy(eigenvectors, mat, ndim*ndim*sizeof(real));
 534     fprintf(stderr, "\nDiagonalizing ...\n");
 535     fflush(stderr);
 536     eigensolver(eigenvectors, ndim, 0, ndim, eigenvalues, mat);
 537     sfree(eigenvectors);
 538
 539     /* now write the output */
 540
 541     sum = 0;
 542     for (i = 0; i < ndim; i++)
 543     {
 544         sum += eigenvalues[i];
 545     }
 546     fprintf(stderr, "\nSum of the eigenvalues: %g (%snm^2)\n",
 547             sum, bM ? "u " : "");
 548     if (fabs(trace-sum) > 0.01*trace)
 549     {
 550         fprintf(stderr, "\nWARNING: eigenvalue sum deviates from the trace of the covariance matrix\n");
 551     }
 552
 553     /* Set 'end', the maximum eigenvector and -value index used for output */
 554     if (end == -1)
 555     {
 556         if (nframes-1 < ndim)
 557         {
 558             end = nframes-1;
 559             fprintf(out, "WARNING: there are fewer frames in your trajectory than there are\n");
 560             fprintf(out, "degrees of freedom in your system. Only generating the first\n");
 561             fprintf(out, "%d out of %d eigenvectors and eigenvalues.\n", end, (int)ndim);
 562         }
 563         else
 564         {
 565             end = ndim;
 566         }
 567     }
 568
 569     fprintf(stderr, "\nWriting eigenvalues to %s\n", eigvalfile);
 570
 571     sprintf(str, "(%snm\\S2\\N)", bM ? "u " : "");
 572     out = xvgropen(eigvalfile,
 573                    "Eigenvalues of the covariance matrix",
 574                    "Eigenvector index", str, oenv);
 575     for (i = 0; (i < end); i++)
 576     {
 577         fprintf (out, "%10d %g\n", (int)i+1, eigenvalues[ndim-1-i]);
 578     }
 579     gmx_ffclose(out);
 580
 581     if (bFit)
 582     {
 583         /* misuse lambda: 0/1 mass weighted analysis no/yes */
 584         if (nfit == natoms)
 585         {
 586             WriteXref = eWXR_YES;
 587             for (i = 0; i < nfit; i++)
 588             {
 589                 copy_rvec(xref[ifit[i]], x[i]);
 590             }
 591         }
 592         else
 593         {
 594             WriteXref = eWXR_NO;
 595         }
 596     }
 597     else
 598     {
 599         /* misuse lambda: -1 for no fit */
 600         WriteXref = eWXR_NOFIT;
 601     }
 602
 603     write_eigenvectors(eigvecfile, natoms, mat, TRUE, 1, end,
 604                        WriteXref, x, bDiffMass1, xproj, bM, eigenvalues);
 605
 606     out = gmx_ffopen(logfile, "w");
 607
 608     time(&now);
 609     gmx_ctime_r(&now, timebuf, STRLEN);
 610     fprintf(out, "Covariance analysis log, written %s\n", timebuf);
 611
 612     fprintf(out, "Program: %s\n", argv[0]);
 613     gmx_getcwd(str, STRLEN);
 614
 615     fprintf(out, "Working directory: %s\n\n", str);
 616
 617     fprintf(out, "Read %d frames from %s (time %g to %g %s)\n", nframes, trxfile,
 618             output_env_conv_time(oenv, tstart), output_env_conv_time(oenv, tend), output_env_get_time_unit(oenv));
 619     if (bFit)
 620     {
 621         fprintf(out, "Read reference structure for fit from %s\n", fitfile);
 622     }
 623     if (ndxfile)
 624     {
 625         fprintf(out, "Read index groups from %s\n", ndxfile);
 626     }
 627     fprintf(out, "\n");
 628
 629     fprintf(out, "Analysis group is '%s' (%d atoms)\n", ananame, natoms);
 630     if (bFit)
 631     {
 632         fprintf(out, "Fit group is '%s' (%d atoms)\n", fitname, nfit);
 633     }
 634     else
 635     {
 636         fprintf(out, "No fit was used\n");
 637     }
 638     fprintf(out, "Analysis is %smass weighted\n", bDiffMass2 ? "" : "non-");
 639     if (bFit)
 640     {
 641         fprintf(out, "Fit is %smass weighted\n", bDiffMass1 ? "" : "non-");
 642     }
 643     fprintf(out, "Diagonalized the %dx%d covariance matrix\n", (int)ndim, (int)ndim);
 644     fprintf(out, "Trace of the covariance matrix before diagonalizing: %g\n",
 645             trace);
 646     fprintf(out, "Trace of the covariance matrix after diagonalizing: %g\n\n",
 647             sum);
 648
 649     fprintf(out, "Wrote %d eigenvalues to %s\n", (int)end, eigvalfile);
 650     if (WriteXref == eWXR_YES)
 651     {
 652         fprintf(out, "Wrote reference structure to %s\n", eigvecfile);
 653     }
 654     fprintf(out, "Wrote average structure to %s and %s\n", averfile, eigvecfile);
 655     fprintf(out, "Wrote eigenvectors %d to %d to %s\n", 1, end, eigvecfile);
 656
 657     gmx_ffclose(out);
 658
 659     fprintf(stderr, "Wrote the log to %s\n", logfile);
 660
 661     return 0;
 662 }