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